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General Notation 


a = perpendicular distance between surfaces for view factor calculation 

b = radius of circular radiating area for view factor calculation 

c « radius of circular absorbing area for view factor calculation 

Cj = coefficients in polynomial or nondimensional prediction functions 

C - thermal equilibrium influence coefficients 

= value of dependent variable (impact damage) at i-th data point 
D = maximum diameter of the bumper hole 
D = minimum diameter of the bumper hole 

MIN 

= average diameter of the hole in the MLI 

D = diameter of projectile 
p 

D pw = average diameter of the pressure wall hole. 

= bumper stand-off distance 

D * estimated value of dependent variable (impact damage) at an 
interpolation or prediction point 
E = elastic modulus of the bumper plate material 

F ■ view factor for radiating from circular area A to circular area A 

A1-A2 12 

F^ ® view factor for radiating from ring area to ring area 

C E b/ a 

G = solar constant 

s 

h = the convective heat transfer coefficient 
h = effective heat transfer coefficient of the dacron netting 

N 

H = c/a 

k = in-plane thermal conductivity of a layer 

L = the distance the module air travels along the pressure wall before 
meeting an obstruction 

M « number of data points in database or number of material properties in 
each record of the materials data file 

N = number of independent variables (impact parameters); or number of 
nodes per layer 

q = the heat flux to the i-th node of the pressure wall 

q^ = net heat flux into the i-th node of a layer 


iv 



q jn = heat flux to a layer at radial position r from adjacent layers 
= heat flux through a layer of dacron netting 
q out ” ^ eat f rom a layer at radial position r to adjacent layers 

q = radiation heat flux 

r 

r = radial position 

2 

R = coefficient of determination 

* distance from i-th data point to interpolation or prediction point 
S = length of influence of a data point 

t = thickness of a layer 

T = temperature in a layer 

T *= bumper thickness 

b 

T « the temperature of the i-th node 

T * temperature of the i-th node of the j-th layer 
hj 

T = pressure wall thickness 

pw 

T = the free stream air temperature of the spacecraft module 
00 

u = the velocity of the module air next to the pressure wall 
00 


V 

V 

1 

X 

X 

J.i 

X 

J.INT 

Ax 

J.l 

A 


= projectile velocity 

/E ' 

= speed of sound in the bumper material = v /p 
« measured value In linear regression equation 
= j-th coordinate (Independent variable) oF i-th data point 
= j-th coordinate (independent variable) of interpolation 
point 

= x - x 

J,l J.INT 

= radial distance between nodes in a layer 


or 


prediction 


c = emissivity of a radiating surface 

<f> = impact angle 

p = mass density of the bumper plate material 

c r = Stef an-Boltzmann constant « 5.6697E-8 W/(m 2 K 4 ) 


0 = weighting factor of a data point 
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Notation for Condensate Modeling 


Cp = specific heat at constant pressure [J/kg K] 

h - enthalpy [J/kg] 

k = thermal conductivity [W/m Kl 

p = pressure [Pa] 

q = heat flux [W/m Z ] 

r « radial coordinate [m] 

rc = radial distance to center of T-cell [m] 

ru = radial distance to center of u-cell [m] 

A r = width of T-cell in radial direction [m] 

= source term for dependent variable <f> 
t = time [s] 

T ® temperature [K] 

u ® velocity in radial direction [m/s] 

Au ■ width of u-cell in radial direction [m] 
v * velocity in axial direction [m/s] 

Av = width of v-cell in axial direction [m] 

z = axial coordinate [m] 

Az = width of T-cell in axial direction 
a = thermal diffusivity [m /s] 

■> coefficient of dynamic viscosity [N s/m ] 
v * coefficient of kinematic viscosity [m Z /s] 

<f> = dependent variable 

p = density [kg/m ] 

T = exchange or diffusion coefficient 

Subscripts 


l .. 

.. radial 

L ... 

. left face of cell 

J ■■ 

.. axial 

R ... 

. right face of cell 

u .. 

.. u-cell 

F ... 

. front face of cell 

V .. 

.. v-cell 

A ... 

. aft face of cell 


Superscripts 

n .... present time level ~ upstream value of variable 
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1. Introduction 

This report describes a computer program called MLITemp that is intended 
to be a design tool for aerospace engineers. The program first uses empirical 
equations to predict hypervelocity impact damage to spacecraft due to space 
debris in earth orbit. A Whipple [1] style of spacecraft wall configuration 
is assumed as is shown in Fig. 1.1. Then, the program predicts the thermal 
effects associated with impact damage, including the amount of condensate 
that would form. 

MLITemp is written in Microsoft BASIC and is designed to run on an 
MS-DOS based personal computer. All of the various capabilities of the 
MLITemp are linked together in the seamless environment of a pull down menu 
system. A help file is provided to assist the user with the menu choices. A 
software user guide is provided in Section 2 of this report. 

Three different techniques for empirically predicting the hypervelocity 
impact damage are provided. An explanation of how each of these techniques 
functions is provided in Sections 3 through 5. More details on these 
empirical prediction techniques are published in a recent NASA Technical 
Memorandum 12). 

The theory behind the thermal analysis program is given in Section 6. 
The thermal analysis methodology that is used in MLITemp was validated by 
experimental testing [3], Also, some thermal system parameters derived during 
the course of this experimental testing are used by MLITemp. 

If the pressure wall of the spacecraft drops below the dew point 
temperature of the spacecraft module air then condensation will tend to 
occur. Such condensate could be hazardous to electrical equipment and could 
also promote the formation of mold. MLITemp estimates the volume of 
condensate that would form. The methodology used to do this is discussed in 
Section 7. 

2. Software User Guide 

2.1 COMPUTER SYSTEM REQUIREMENTS 

The software developed for this project was written using the Microsoft 
BASIC Professional Development System (BPDS). However, the programs of 
MLITemp that do not use the menu, window, and mouse toolbox of BPDS can be 
modified and recompiled using Microsoft QuickBASIC or some other language. An 
EGA or VGA graphics card and monitor, and an Intel 80286, 80386 or 80486 CPU 
is required to run the software. A math coprocessor must be available. The 
software is provided on 360K computer disks. 
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2.2 PROGRAM AND DATA FILES OF MLITEMP 

An annotated listing of the program and data files of MLITemp follows: 


MLITEMP. BAS - source code for the main program that runs the other programs 
(ASCII). 

MLITEMP.EXE - compiled version of the main program. 

DATABASE. BAS - source code for the database creation program (ASCII). 

DATABASE.EXE - compiled version of the database creation program. 

DBASEDEL. BAS - source code for the database record deletion program (ASCII). 

DBASEDEL.EXE - compiled version of the database record deletion program. 

DBASEOUT.BAS - source code for the database viewing program (ASCII). 

DBASEOUT.EXE - compiled version of the database viewing program. 

EDITIMPA.BAS - source code for the impact parameters editing program 

which operates on the impact parameters file impact, par 

(ASCII). 

EDITIMPA.EXE - compiled version of the impact parameters editing program. 

EDITTHER.BAS - source code for the thermal parameters editing program 

which operates on the thermal parameters file thermal. par 

(ASCII). 

EDITTHER.EXE - compiled version of the thermal parameters editing program. 

INVRMETM.BAS - source code for the inverse R method damage prediction program 
(ASCII). 

INVRMETH.EXE - compiled version of the inverse R method damage prediction 
program. 

POLYMETH.BAS - source code for the polynomial function damage prediction 
program (ASCII). 

POLYMETH.EXE - compiled version of the polynomial function damage prediction 
program. 

NONDIMEN. BAS - source code for the nondimensional function damage prediction 
program (ASCII). 

NONDIMEN.EXE - compiled version of the nondimensional function damage 
prediction program (ASCII). 

PREVNOND.BAS - source code of the program that prompts the user if 
previously calculated nondimensional function coefficients 
should be used in the calculations (ASCII). 

PREVNOND.EXE - compiled version of the prevnond.bas program. 

SHOWIMPA. BAS - source code of the program that displays impact predictions on 
the screen (ASCII). 

SHOWIMPA.EXE - compiled version of the showlmpa.bas program. 

UPDATE. BAS - source code of the program that updates the thermal parameters 
file with the latest impact results data (ASCII). 

UPDATE.EXE - compiled version of the update. bas program. 

THERMAL. BAS - source code of the program that performs the thermal analysis 
(ASCII). 

THERMAL. EXE - compiled version of the thermal. bas program. 

CONDEN.BAS - source code of the program that performs the condensation 
calculations (ASCII). 

CONDEN.EXE - compiled version of the program that performs the condensation 
calculations. 

SHOWTEMP. BAS - source code of the program that displays the results of the 
thermal analysis on the computer screen. 

SHOWTEMP.EXE - compiled version of showtemp.bas. 

HELP. BAS - source code of the program that displays the help file help.doc. 

HELP.EXE - compiled version of the help. bas program. 

MATERIAL.DAT - a typical database file of material properties which is used 
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by the 1NVRMETH program (ASCII). 

MLI.DAT - a typical database file of experimental results (ASCII). 

IMPACT.PAR - a typical impact parameters file (ASCII). 

THERMAL. PAR - a typical thermal parameters file (ASCII). 

NONDIMEN. OLD - a file storing the previously calculated nondimensional 
function coefficients (ASCII). 

UPDATE. PAR - a file storing the coefficients of the diameter ratio function 
used by the update.exe program (ASCII). 

CONDEN.PAR - a data file used for transferring information from the 
thermal.exe program to the conden.exe program (ASCII). 

HELP.DOC - this Is the help document displayed by the help.exe program. 

2.3 SOFTWARE INSTALLATION AND EXECUTION 

The software is installed by first creating a subdirectory on the hard 
disk and then copying all of the files from the computer disks into that 
subdirectory. If disk space is a problem then the source code files 
(/tlename.BAS.) need not be copied. The program is started by typing MLITEMP. 
The options of MLITEMP can be selected from the keyboard or by using the 
mouse as will now be described. 

r 

WARNING - Be sure you are using the correct units / The correct units for the 
various data files are given In Section 2.4. 

The standard procedure for running MLITemp Is as follows. First, the 

impact parameters file (IMPACT.PAR) is edited to reflect the desired 
hypervelocity impact conditions. Then, the hypervelocity impact testing 
results file Is edited if necessary. A typical impact testing results file 
called MLI.DAT is provided on disk. Next, one of the three impact damage 

prediction programs (inverse R, polynomial fit, nondimensional function) are 
run. The thermal parameters file (THERMAL.PAR) is then updated with the 

impact damage prediction results and possibly edited with respect to other 
thermal properties. Finally, the thermal analysis program Is executed and the 
results viewed. More details on these procedures are provided below. The 
procedures described above can be performed by selecting tasks from the menu. 
The menu can be activated by clicking with the mouse, or by pressing the 

<ALT> key. Menu commands can be selected by using the mouse, by using the 
arrow keys and pressing <ENTER>, or by typing the red letter of each command. 
The menu commands of MLITemp are described below. 

Main Menu - FILE 

ADD TO IMPACT DATA FILE: 

This allows the user to add data records to the hypervelocity impact 
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testing results file. This file contains the data used for making empirical 
predictions of impact damage to the bumper, the multilayer insulation (MLI), 
and to the pressure wall. This menu pick runs the program DATABASE.EXE. The 
first record in the impact parameters file, IMPACT.PAR, contains the name of 
the impact testing results file that will be operated on. This filename can 
be changed by selecting the EDIT IMPACT PARAMETERS FILE menu pick, which is 
described below. The user can move from edit box to edit box in the editor 
window by pressing the <TAB> key, by pressing the <ENTER> key, or by pressing 
the up or down arrow keys. The user can move around within an edit box of an 
edit window using the <HOME>, <END>, and arrow keys. A button at the bottom 
of the edit window (add to database, cancel this data entry, exit program) 
can be activated by pressing <ENTER> after a button has been selected using 
the <TAB> or arrow keys. The mouse can also be used to move between edit 
boxes and buttons. A data input window is shown In Fig. 2.3.1. Note that the 
program automatically inserts defaults for data values that seldom vary. 

REMOVE FROM IMPACT DATA FILE: 

This allows the user to remove records from the hypervelocity impact 
testing results file. This file contains the data used for making empirical 
predictions of impact damage to the bumper, the MLI, and to the pressure 

wall. This menu pick runs the program DBASEDEL.EXE. The first record in the 
impact parameters file, IMPACT.PAR, contains the name of the Impact testing 
results file that will be operated on. This filename can be changed by 

selecting the EDIT IMPACT PARAMETERS FILE menu pick, which is described 

below. A button at the bottom of the edit window (OK to remove, quit) can be 
activated by pressing <ENTER> after a button has been selected using the 
<TAB> or arrow keys. The mouse can also be used to move between buttons. The 
data window of DBASEDEL.EXE is shown in Fig. 2.3.2. 

VIEW IMPACT DATA FILE: 

This allows the user to view records In the hypervelocity impact testing 
results file. This file contains the data used for making empirical 
predictions of impact damage to the bumper, the MLI, and to the pressure 

wall. This menu pick runs the program DBASEOUT.EXE. The first record in the 
impact parameters file, IMPACT.PAR, contains the name of the impact testing 
results file that will be operated on. This filename can be changed by 

selecting the EDIT IMPACT PARAMETERS FILE menu pick, which is described 

below. A button at the bottom of the edit window (next data record, exit 
program) can be activated by pressing <ENTER> after a button has been 

selected using the <TAB> or arrow keys. The mouse can also be used to move 
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between buttons. A typical data record view window is shown in Fig. 2.3.3. 

EDIT IMPACT PARAMETERS FILE: 

This allows the user to view and edit the impact parameters file, 

IMPACT. PAR. This menu pick runs the program EDITIMPA.EXE. The user can move 
from edit box to edit box in the editor window by pressing the <TAB> key, by 

pressing the <ENTER> key, or by pressing the up or down arrow keys. The user 

can move around within an edit box of an edit window using the <HOME>, <END>, 
and arrow keys. A button at the bottom of the edit window (save changes and 
exit, exit program) can be activated by pressing <ENTER> after a button has 
been selected using the <TAB> or arrow keys. The mouse can also be used to 
move between edit boxes and buttons. The impact parameters edit window is 
shown in Fig. 2.3.4. 

EDIT THERMAL PARAMETERS FILE: 

This allows the user to view and edit the thermal parameters file, 

THERMAL. PAR. This menu pick runs the program EDITTHER.EXE. The user can move 
from edit box to edit box in the editor window by pressing the <TAB> key, by 

pressing the <ENTER> key, or by pressing the up or down arrow keys. The user 

can move around within an edit box of an edit window using the <HOME>, <END>, 
and arrow keys. A button at the bottom of the edit window (next window, save 
changes and exit, exit program) can be activated by pressing <ENTER> after a 
button has been selected using the <TAB> or arrow keys. Four windows are 
required to view all of the thermal parameters. The user can proceed from 
window to window using the <NEXT WINDOW> button. The mouse cafi also be used 
to move between edit boxes and buttons. Fig. 2.3.5 Illustrates the four 

windows of the thermal parameters editing program. 

CURRENT DIRECTORY FILENAMES: 

This menu pick causes the names of the files in the current directory to 
be listed on the screen. This may be useful if the user forgets the name of a 
data file. 

DOS SHELL: 

This menu causes a DOS shell to be created. This will allow the user to 
copy files and perform other tasks without leaving the MLITemp program 
permanently. Entering "exit" causes the DOS shell to close. 

EXIT: 

This menu pick will end the MLITemp program. 

Main Menu - IMPACT 

INVERSE R METHOD: 

This menu pick will cause a hypervelocity Impact damage prediction to be 
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made using the "inverse - R" prediction algorithm (program INVRMETH.EXE is 
executed). The details of this prediction algorithm are described in Section 
3. The empirical prediction is based on experimental data contained in the 
impact testing results file named in the impact parameters file (IMPACT. PAR). 
The impact parameters associated with the prediction are contained in file 
IMPACT. PAR. File IMPACT. PAR can be edited from the EDIT IMPACT PARAMETERS 
FILE menu pick under the FILE main menu. The bumper hole major and minor 
diameters, the MLI hole diameter, and the pressure wall hole diameter are 
predicted. 

POLYNOMIAL FIT METHOD: 

This menu pick will cause a hypervelocity impact damage prediction to be 
made using the "polynomial fit" prediction algorithm (program POLYMETH.EXE is 
executed). This prediction algorithm is described in Section 4. The empirical 
prediction is based on experimental data contained in the impact testing 
results file named in the Impact parameters file (IMPACT.PAR). Hie impact 
parameters associated with the prediction are contained in file IMPACT.PAR. 
File IMPACT.PAR can be edited from the EDIT IMPACT PARAMETERS FILE menu pick 
under the FILE main menu. The bumper hole major and minor diameters, the MLI 
hole diameter, and the pressure wall hole diameter are predicted. 
NONDIMENSIONAL FUNCTION METHOD: 

This menu pick will cause a hypervelocity impact damage prediction to be 
made using the "nondimensional function" prediction algorithm (program 
NONDIMEN.EXE is executed). Details of this prediction algorithm are given in 
Section 5. The empirical prediction is based on experimental data contained 
in the impact testing results file named in the impact parameters file 
(IMPACT.PAR). The impact parameters associated with the prediction are 
contained in file IMPACT.PAR. File IMPACT.PAR can be edited from the EDIT 
IMPACT PARAMETERS FILE menu pick under the FILE main menu. The bumper hole 
major and minor diameters, the MLI hole diameter, and the pressure wall hole 
diameter are predicted. This program takes a relatively long time to run 
since 23 nonlinear function coefficients are being fit to the experimental 
data. At the end of the execution these coefficients are stored in a file 
named NONDIMEN. OLD. The following scheme was developed to speed up the 
calculations for the case where there has been no change in the Impact 
testing results file (and thus the function coefficients would not change). 
Before running the NONDIMEN.EXE program a program called PREVNOND.EXE is run. 
PREVNOND.EXE prompts the user for whether the old nondimensional function 
coefficients should be used. If the user picks yes then the contents of file 
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NONDIMEN. OLD are copied to a file called NONDIMEN.NEW. If the program 
NONDIMEN.EXE senses the existence of file NONDIMEN.NEW, then no new function 
coefficients are calculated, and the old function coefficients (that were 
originally contained in file NONDIMEN. OLD) are used to make the damage 
predictions. 

SHOW CURRENT IMPACT RESULTS: 

This menu pick displays the current set of Impact predictions (program 
SHOWIMPA.EXE is run). The Impact predictions are stored in a file named in 
the impact parameters file IMPACT. PAR. Also, SHOWIMPA.EXE is automatically 
run after each of the damage prediction programs have completed their 
calculations. Typical output from this program is shown in Fig. 2.3.6. Press 
<ENTER> to exit from this program. 

UPDATE THERMAL PARAMETERS FILE: 

This menu pick (which executes program UPDATE.EXE) updates the thermal 
parameters file, THERMAL.PAR, with the bumper and MLI hole diameters obtained 
from the most recent run of an impact damage prediction program. UPDATE.EXE 
performs three operations. First, it determines an average bumper hole 
diameter by averaging the major and minor bumper hole diameters calculated by 
an impact prediction program. Then, UPDATE.EXE converts the average diameter 
from units of Inches to units of meters, as required by the thermal analysis 
program. Finally, at the option of the user, the MLI hole diameter is 
adjusted with the "diameter ratio" parameter. The diameter ratio parameter Is 
an empirical function of the Impact parameters. The diameter ratio parameter 
is intended to account for the fact that the apparent MLI hole diameter from 
a thermal analysis context is in general different from the observed MLI hole 
diameter. The six empirical function coefficients used to calculate the 
diameter ratio parameter are stored in ASCII file UPDATE.PAR which may be 
modified by the user. Details on the diameter ratio function are given in 
Section 6. 

Main Menu - TEMPERATURE 

PERFORM THERMAL CALCULATIONS: 

This menu pick runs the thermal analysis program THERMAL.EXE. The theory 
behind the thermal calculations is described in Section 6. The analysis is 
based on parameters contained In the thermal parameters file, THERMAL.PAR. 
The initial values used f or the thermal analysis and the results of the 
thermal analysis are sent to files named In THERMAL.PAR. File THERMAL.PAR can 
be edited f rom the EDIT THERMAL PARAMETERS FILE menu pick under the FILE main 
menu. If pressure wall temperatures drop below the dew point temperature of 
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the spacecraft module air, then the condensate thickness profile is 

calculated by program CONDEN.EXE. Details on the condensate calculation 
program are given in Section 7. The thermal analysis program, THERMAL.EXE, 
and the condensate analysis program, CONDEN.EXE, compiunicate to each other by 
means of the data file CONDEN.PAR. After the thermal and condensate analysis 
has been completed, the results are graphically illustrated on the screen by 
program SHOWTEMP.EXE. The results are also written to a file named in 
THERMAL. PAR. A typical thermal results file is shown in Table 2.3.1. 

SHOW CURRENT THERMAL RESULTS: 

This menu pick runs program SHOWTEMP.EXE which illustrates the results 
of thermal and condensate analyses on the computer screen. Color contour 
plots of the bumper and pressure wall temperature distributions are shown. 
Also, a cross section through the geometric configuration of the modeled 
section of the spacecraft wall and the condensate layer (if present) are 
drawn to scale on the screen. A typical display of results is illustrated in 
Fig. 2.3.7 (color contours can not be seen in the figure). 

Main Menu - HELP 

VIEW HELP DOCUMENT: 

This menu pick will cause program HELP.EXE to run which displays an 
ASCII file containing instructions on how to use MLITemp. Additional 
information may be added to this file by using a text editing program if the 
user so desires. 

2.4 DESCRIPTION OF THE DATA FILES OF MLITEMP 
File - MATERIAL.DAT 

The MATERIAL.DAT file that is provided on disk as an example of a 
typical materials data file. Any valid DOS name can be used for this file. 
Thus, the user may have several of this type of data file in a directory for 
different purposes. A file of this nature Is required while running the 
inverse R program. The materials data file is an ASCII file that can be 
created and modified using any standard text editor. The format of the file 
is as follows: 


■ material property 1 name string 
m material property 2 name string 


LISTING OF NAMES OF MATERIAL 
► PROPERTIES TO BE MODELED (MAXIMUM OF 10) 
(25 CHARACTERS MAX) 


■ material property M name string 

( 

m material 1 name string 
m material property 1 for material I 
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■ material property 2 for material 1 


► TYPICAL DATA RECORD 


■ material property M for material 1 

, ► ANY NUMBER OF DATA RECORDS MAY BE USED 


A material data file provided with the MLITemp software is called 

MATERIAL.DAT and is reproduced below: 

Density (lb/in~3) 

Elastic Mod , (lb/irC2) 

Ultimate Strgth (lb/in"2) 

Sp . Heat (BTU/(lb-deg R)) 

Melting Temp (deg R) 

{ 

1100 

9.780E-2 

1.000E7 

1.600E4 

2.140E-1 

1.680E3 

) 

( 

2219-T87 

1.030E-1 

1.050E7 

6.300E4 

2.050E-1 

1.680E3 

} 

( 

6061-T6 
9.800E-2 
9.900E6 ' 

4.200E4 

2.100E-1 

1.680E3 

) 

The MATERIAL.DAT file listed above is set up to model the material 
properties: density, elastic modulus, ultimate strength, specific heat, and 

melting temperature. Other physical properties can be used to a maximum of 
10. The units do not have to be included in the material property name 
string. MATERIAL.DAT contains three records of material data for materials: 
1100, 2219-T87, and 6061-T6. The material names are treated as string 
variables and thus can be any combination of numbers and letters. Any number 
of records of material data may be included. The order of the material 
properties must be the same In every record and must be ordered as the 
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material property name strings are listed. For instance, referring to file 

MATERIAL.DAT, the specific heat of material 2219-T87 is 2.050E-1. 

The purpose of the material properties database file is to provide an 
efficient, yet very flexible scheme for inputting material property data into 
the inverse R method computer program. The user can easily change the 
material properties to be modeled without disturbing the database file of 
experimental results. If the materials used for the projectile, bumper and 

pressure wall do not vary in the database, then the contents of the material 
properties database file will have no effect on the damage predicted by the 
Inverse R method program. The polynomial function method program assumes that 
the material properties do not vary in the database. The nondimensional 
function method program assumes that the material properties of the 
projectile and pressure wall do not vary in the database and assumes the 

bumper Is constructed of an aluminum alloy. 

HYPERVELOCITY IMPACT TESTING RESULTS FILE 

The other database file required for running the Impact damage 
prediction programs is associated with the experimental data. This file can 
be created (and enlarged) by running the database maintenance programs 
described in the previous section or it can be created using any standard 

text editor since it is an ASCII file. This file can be given any valid DOS 
file name. Currently, up to 100 data records can be placed, in this file. The 
format for this file is as follows: 

{ 

■ test ID number 

a source of the data 
a test date 

m bumper material name - which must be of the same format as that listed in 
the material data file 
a bumper thickness (inches) 
a bumper standoff (Inches) 

m pressure wall material name - which must be of the same format as that 
listed in the material data file 
a pressure wall thickness (Inches) 

■ projectile material name - which must be of the same format as that listed 
in the material data file 

m projectile diameter (inches) 

m impact angle (degrees) - this is the angle between the normal to the bumper 
and the line of travel of the projectile 
a projectile velocity Ckm/sec) 
m major axis of bumper hole (Inches) 
m minor axis of bumper hole (inches) 
m average Mil hole diameter (inches) 
a average pressure wall hole diameter (Inches) 

) 
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MLI.DAT Is provided as an example of an experimental database file. This file 
is stored on the computer disks. It contains information on the specimens 
recently used for thermal testing in Sunspot Thermal Vacuum Chamber of MSFC. 
To help understand the format information given above, the first record of 
MLI.DAT is presented below for comparison: 

( 

1012 

MSFC 

05/08/90 

6061-T6 

,08 

4 

2219-T87 

.125 

1100 

.313 

0 

6.72 

.729 

.729 

2.2 

.375 

File - IMPACT.PAR 

IMPACT. PAR specifies the following information related to predicting the 
damage associated with an impact event (required units shown in brackets): 


■ filename of the hypervelocity impact testing results file 
m filename of the material data file 

m filename of the file to be used to store the output from the impact damage 
prediction programs - the impact damage prediction programs can be accessed 
under the IMPACT main menu 

■ bumper material name - which must be of the same format as that listed in 
the material data file 

m bumper thickness (inches) 

■ bumper standoff (inches) 

m pressure wall material name - which must be of the same format as that 
listed in the material data file 
m pressure wall thickness (inches) 

m projectile material name - which must be of the same format as that listed 
in the material data file 
m projectile diameter (inches) 

m impact angle (degrees) - this is the angle between the normal to the bumper 
and the line of travel of the projectile 

■ projectile velocity (km/sec) 


File - THERMAL. PAR 

THERMAL. PAR specifies the following information related to predicting 
the temperature in the spacecraft walls and also condensate layer 


11 



thickness (required units shown in brackets); 


■ filename of the file to be used to store the output from the thermal 
analysis program - the thermal analysis program is under the TEMPERATURE 
main menu 

m filename of the file uses to store initial temperatures and heat fluxes for 
the thermal analysis model - if the Initial values file does not exist then 
estimated pressure wall and bumper temperatures will be used 
m MLI hole diameter (m) 

m MLI standoff (m) - this is the distance from the outer surface of the 
pressure wall to the centerline of the MLI blanket 
a estimated pressure wall temperature (K) - this is only used if the initial 
values file does not exist 

a estimated bumper temperature (K) - this is only used if the initial values 
file does not exist 

a temperature conversion factor 1 - used for converting from degrees K to 
degrees F for display purposes 

a temperature conversion factor 2 - used for converting from degrees K to 
degrees F for display purposes 

■ number of aluminized MLI layers 

a radius of the area thermally modeled (m) - uniform temperatures are assumed 
to exist outside of this area 
a pressure wall thickness (m) 
a thickness of an aluminized MLI layer (m) 
a beta cloth thickness (m) 
a bumper thickness (m) 
a bumper standoff (m) 

a space thermal radiation - influx from far- field radiator (W/m~2) - if the 
area of interest is facing In the direction of deep space then this 
parameter should have a magnitude of zero - if the area of interest is 
facing the sun then this parameter should have a magnitude of 1353 W/m*2 
(known as "solar constant ", Gs) 
a pressure wall thermal conductivity (W/mK) 

a MLI aluminized layer thermal conductivity (W/mK)* dacron netting heat 
transfer coefficient (W/m~2K) - recent experiments have shown that this 
parameter should have a value of 1.0687 W/m*2K for baselined Space Station 
MLI 

a beta cloth thermal conductivity (W/mK) 
a bumper thermal conductivity (W/mK) 
a pressure wall emissivity 
a aluminized MLI layer emissivity 

a beta cloth emissivity : =- 

■ bumper emissivity outward - this allows for a special coating on the 
outside of the bumper 

a bumper emissivity inward 

a Stef an-Boltzmann constant = 5.6697E-8 W/m"2K"4 
a maximum number of global Iterations for a given mesh size 
a convergence factor 

a initial number of nodes in each layer of the model 
m final number of nodes in each layer of the model 

■ bumper hole diameter (m) 

m mean air temperature of the spacecraft module Interior (K) 

■ spacecraft module air dew point temperature (K) - typically 290 K 
m spacecraft wall convective heat transfer coefficient (W/m*2K) 

m condensate density (kg/m"3) 
a condensate kinematic viscosity (m*2/s) 
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■ condensate thermal conductivity (W/mK) 

m condensate constant pressure specific heat (J/kgK) 
m spacecraft module air density ( kg/m~3 ) 

■ spacecraft module air kinematic viscosity (m~2/s) 
m spacecraft air thermal conductivity (W/mK) 

m spacecraft air constant pressure specific heat (J/kgK) 

In the next three sections, the three techniques used for predicting 
Impact damage are discussed. 

3. The Inverse R Prediction Technique 

The usual procedure for making predictions from experimental data is to 
assume some form for the equation relating the independent variables to the 
dependent variable. A function of this nature is described in Section 5 of 
this report. The equation typically contains empirical coefficients, the 

values of which are determined from a fit to the experimental data [4-9], The 
method of least squares (maximizing the coefficient of determination, R"2) is 
an example of a popular technique for obtaining the coefficients from the 
experimental data. The final result is a closed form equation for making 
predictions. 

This approach has been found to work very well for many engineering 
applications, however there are some disadvantages. A suitable form for the 
prediction equation must be developed. This is often difficult. Incorporating 
additional independent variables in an existing equation can pose problems. 
Usually, a well defined procedure for taking Into account new experimental 
data is not put in place. Generally, a single set of empirical coefficients 

are used to make predictions over a fairly wide range of values of the 

independent variables. Thus, the best data in a database for making a 

prediction with a particular set of independent variables may not be used to 
best advantage. Also, it is usually difficult to assess the accuracy of a 
particular prediction. 

In this section, a new method (called inverse R method) for making 

empirical predictions based on experimental data Is discussed. The method 
uses a very general form of prediction equation that can be applied in the 

same manner to all problems. Thus, the user is not required to develop a 
suitable form for the prediction equation and additional independent 

variables can be easily incorporated. The new method is designed to work off 
a database that can be continuously updated as new experimental data becomes 
available. The method automatically takes advantage of the most appropriate 
data in the database for a given set of independent variables. 

The new technique consists of four main steps which will now be 
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described. 

Step 1. Normalize the Independent Variables. 

In general, the independent variables will vary greatly in magnitude. In 
hypervelocity impact work, dimensions can be of order 10 and velocities of 
order 10 6 . The new technique requires that all variables be of the same order 
of magnitude. This was accomplished by scaling the independent variables such 
that their mean value was equal to unity. Other scaling methods could perhaps 
be used to improve the accuracy of this technique. For instance, the 
variables could be scaled such that predicted values of points in the 
database more closely match the measured values. This scaling technique was 
not tested. The dependent variables need not be scaled. 

This technique works off a database that can and should be kept updated 
with the latest experimental data. Thus, the scaling factors will change as 
time progresses and the size of the database increases. 

Step 2. Select a Series of Points In the Data Domain For Interpolation. 

Two general requirements for prediction schemes are: the method should 
be capable of smoothing the data to (hopefully) cancel out the random scatter 
typically present in experimental measurements, and the technique should 
allow for making reliable predictions outside of the domain of the measured 

data. Here, these requirements are satisfied by using the data to make ten 

Interpolations from within the domain of the data, which are then used for 
predicting the dependent variable at some point of interest. The ten 
"interpolation" points should provide for sufficient smoothing of the data 
and also capture the trend characteristics of the data for extrapolation 

purposes, if an extrapolation is required. The number of interpolation points 
to use was selected on the basis of trial and error. Note, in some cases 
extrapolation can produce misleading results regardless of the extrapolation 

technique used. 

Fig. 3.1 provides an illustration of how the interpolation points are 
selected for a hypothetical case with two independent variables. An identical 
approach is used for the case of an arbitrary number of independent 
variables. In Fig. 3.1, the independent variables are in the plane of the 
page and the dependent variable takes the form of a surface out of the plane 
of the page. 

First, a prediction "vector" is drawn from the origin through the point 
in the domain where a prediction of the dependent variable Is required, which 
is called the "target" point. Then the "min" and "max" points (Fig. 3.1) are 
located on the prediction vector by considering the intersection points of 
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perpendiculars from the data points to the prediction vector. The closest 
intersection point to the origin defines the min point, and that of the 
farthest, the max point. Ten equally spaced points (interpolation points) on 
the prediction vector between the min and max point are then used for the 
next step in the prediction process. If the target point lies between the min 
and max points then an interpolation is required, otherwise an extrapolation 
is required. 

Step 3. Estimate Values of the Dependent Variable at Interpolation Points, 

Next, values for the dependent variable must be estimated at the ten 
interpolation points. This is done as indicated in the following equation: 



The distances, R J# are determined by the usual formula for determining the 
"distance” between two points in an N dimensional space: 



where x and x are the J-th coordinates (bumper thickness and so on) 
j,i j*int 

of the data point and the point to be predicted, respectively. The need for 
scaling the independent variables is evident from considering the form of Eq. 
(3.2). 

The form of Eq. (3.1) will now be considered. It is assumed that if all 
measured data points are the same "distance" R from an interpolation point 
then ail the measured data should be given equal weight. This situation is 
illustrated for the case of two independent variables (N * 2) in Fig. 3.2. 
This can be interpreted as saying that each data point has some 

"characteristic length of influence", S, that subtends an angle 0 «* S/R « 
S/R N ~ l as indicated in Fig. 3.2. The 0 can be taken as the weighting factor. 
For the constant R case shown in Fig. 3.2, all data points would be given the 
same weight. Fig. 3.3 illustrates the case for which the data points are 
considered to be equally valid (same S), but are located different 
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"distances" from the Interpolation point. Here, the weighting factors will be 
of the form 6 ( = S/R 1 ^ \ and thus data points closer to the interpolation 
point will be given a higher weight. The value of the dependent variable at 
the interpolation point can be estimated from D = EB^/ZG^ which leads to 
Eq. (3.1) and hence this technique is given the name inverse R method. Note 
that a value for S is not required as it cancels out of the equation. 

The three dimensional (three independent variables) application of this 
procedure leads to equations identical in form to those used for determining 
view factors in the field of radiation heat transfer [10]. Thus, the 
rationale for the inverse R method can be interpreted as follows. The 
measured data points are "radiating" information to the interpolation point. 
The farther the data point is away, the weaker the "radiation” (lower weight 
given to the information). In principle, the method can easily be extended to 
any number of Independent variables, N. 

Step 4. Fit a Polynomial Through the Interpolation Points and Make 
Prediction. 

The final step in the process involves fitting a polynomial through the 
ten interpolation points and then using the polynomial to make a prediction 
of the dependent variable at the target point. The polynomial describes how 
the dependent variable behaves as a function of distance along the prediction 

vector. By trial and error it was found that a fourth-order polynomial worked 
well for this application. The polynomial could be used for interpolation or 

extrapolation depending on the location of the target point. There would of 
course be considerably more uncertainty in the prediction for the case of 
extrapolation. Errors in the ten interpolation points tend to get smoothed by 

the polynomial. 

Because of it’s uniqueness, the inverse R method was tested to ensure it 
would provide reliable predictions. These tests are reported in [2, 11], 

4. The Polynomial Function Prediction Technique 

In this section the polynomial function prediction technique is 

described. This method is based on the concepts associated with the finite 
element method (FEM). In FEM, relatively low order polynomials are used to 

interpolate the functions of interest (such as displacements, temperatures, 
and velocities) over a small portion of domain where the function is active 
called an element. The coefficients of the polynomial are derived from known 

values of the function of interest at points called nodes on the boundary of 

the element. For this application, the nodal values of the functions of 
interest (bumper hole size and so forth) were measured experimentally and are 
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thus known quantities. This technique involves selecting a sufficient number 
of experimental data (node) points and then determining the coefficients of 
the polynomial from this data. 

Ideally, the nodes "closest" to the prediction point in impact parameter 
space should be used to evaluate the polynomial coefficients and thus make a 
prediction. However, the set of closest nodes may not form linearly 

Independent set of data, making it impossible to solve for the polynomial 
coefficients. Thus, remoter nodes must be considered in an attempt to find a 
linearly independent set of data. The technique used for selecting remoter 
nodes is discussed below. In general, the impact parameters will vary greatly 
in magnitude. In hypervelocity impact work, dimensions can be of order 10 and 
velocities of order 10 6 , This polynomial function approach requires a 

reasonable scheme for determining "distances" between data points in impact 
parameter space. This is accomplished in the program by scaling the impact 
parameters (bumper thicknesses and so on) such that their mean value is equal 
to unity. Of course, the dependent variables, such as bumper hole size, need 
not be scaled. Having scaled the independent variables, the usual formula for 
determining the "distance", Rj, between two points (experimental data point 
and the prediction or interpolation point) in a multidimensional space can be 
used: 



1 


N 



J=1 



(4.1) 


where x and x are the J-th coordinates (bumper thickness and so on) 

J,i J»int 

of the data point and the point to be predicted, respectively. The need for 
scaling the independent variables is evident from considering the form of Eq. 
(4.1). 

The form of the polynomial will now be considered. FEM theory dictates 

that a "complete" polynomial should produce the best results [12]. Here we 

have six independent variables (bumper thickness and so forth), x (j = 1 

J* * 

to 6), associated with the i-th experimental data point to consider. It was 

decided to use Ax (= x - x ) values In the polynomial equation to 

J.i J.i J.int 

simplify the calculations. The lowest order complete polynomial for this case 
Is: 


D = C + C *Ax + C *Ax + C «Ax + C *Ax + C *Ax, + C *Ax, (4.2) 

1 1 2 1,1 3 2,1 4 3,1 5 4,1 6 5,1 7 6.1 
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Seven linearly independent data points, D., are required to determine the 
seven polynomial coefficients, C. Eq. (4.2) allows for a linear variation in 
damage along each coordinate axis in the design space. Obviously, allowing 
for a quadratic variation in the damage would provide a much better fit to 
the data. Unfortunately, a ’’complete" quadratic function with six variables 
would require too many linearly independent experimental data points to be of 
practical use considering the relatively small quantity of experimental data 
available. 

Coefficient is the prediction of the damage at the point in the 
design space where the prediction is required, since this is the value jof the 
polynomial (Eq. 4.2) when all Ax are set equal to zero. If one or more of 
the prediction parameters, such as bumper thickness, does not vary in the 
experimental database file then program POLYMETH will sense this and 
automatically take that variable or variables out of Eq. (4.2). If one impact 
parameter does not vary, only six polynomial coefficients need be determined 
and thus only six linearly independent data points are required. 

The method used to select the linearly independent set of data points 
from the database for determination of the function coefficients, C^, of Eq. 
(4.2) will now be discussed. For illustration purposes, assume that three 
independent variables are active and thus four linearly independent data 
points are required to fit coefficients through C 4 * First, the four 

closest data points are selected and tested for linear independence. If they 
are linearly independent, then the coefficients can be determined and the 
prediction made. If the four closest data points are not linearly 
independent, then groups of four data points (the closest data point plus 
three others) are selected from the closest five data points and tested for 
linear independence. The first linearly Independent set of data points found 
is used for coefficient determination. If a set of suitable data points Is 
not found, then sets of four data points are selected from the closest six 
data points and so on. 

The number of ways to choose r items from n items, C(n,r), is given by 
the following equation: 


n! 

C(n, r) ■ (4.3) 

(n-r)!r! 

From Eq. 4.3, there are 20 ways to choose 3 items from 6 items. Thus, as 
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shown in Table 4.1, twenty sets of data would have to be tested for linear 
independence when selecting four point data sets (the closest plus three 
other data points) from the closest 7 data points. Note in Table 4.1 that the 
closest data sets are tested first and data point 1 is always used. 

The effectiveness of this prediction technique is tested in [2]. 

5. The Nondimensional Parameter Prediction Technique 

In many applications It has been found that empirical functions are best 
represented in terms of nondimensional parameters. Reynolds number is an 
example of a nondimensional parameter that has found widespread use in 
empirical equations of fluid mechanics. Program NONDIMEN uses a series of 
empirical functions based on nondimensional parameters of the form given in 
113): 
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The function coefficients were determined using an optimization routine 

to adjust the values of the coefficients so as to maximize the coefficient of 
2 

determination (R ) of each of the functions. Thus, the nondimensional 

functions were adjusted to match the experimental results as closely as 
possible in a least squares sense. This approach to coefficient evaluation is 
suitable for any form of prediction function - linear or nonlinear. The 

nature of the optimization routine will now be described. 
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The magnitudes of the function coefficients can vary by several orders 
of magnitude. To avoid numerical problems it is advisable to work with 
percentage changes in the function coefficients. This approach also provides 
a simple way of controlling the amount of change in the function coefficients 
from one optimization iteration to the next. If the maximum allowable 

percentage change is too large, the optimizer could thrash back and forth 
around the optimum design point without ever converging to it. Alternatively, 

if the maximum allowable percentage change is too small, then it could take 
an impractical number of iterations to get to the optimum design point, or 

the optimizer could get "stuck" in a local maximum of the coefficient of 

determination function before getting to the global maximum. 

The maximum allowable percentage change in the nondimensional function 
coefficient magnitudes used is 1.0 (equivalent to a 1007, change). The 

optimizer is designed to reduce the magnitude of the search domain parameter 
as the optimization process proceeds. The final value will be 1/100 of the 
initial value. The idea here is to allow large changes in the design 
variables initially, to quickly get into the vicinity of the global maximum 
in the design space, and then use finer steps to precisely locate the global 
maximum. The user is free to change this parameter to attempt to improve 
optimization efficiency. 

The initial values of the function coefficients are set equal to zero. 

Optimal values of the function coefficients could be positive* negative or 

zero. 

The method chosen here for search vector selection is based on Powell's 
method 1 141. This Is a first order method that does not require the 
calculation of the gradient vector. Here, Powell’s method was modified as 

follows. Initially, a number of search vectors equal to the number of 
function coefficients are created. The components of these vectors are random 
numbers between -1 and +1. The components of each random search vector are 
then scaled, such that the largest component has a magnitude of unity. These 
vectors are stored as columns of a "search matrix". Next, the coefficient of 
determination is evaluated at the current point in the design space and at 
design points given by ♦/- the search domain parameter times the first column 
of the search matrix. If either of the + or - design points has a coefficient 
of determination greater than that of the current design point, then the 
design point corresponding to the highest coefficient of determination will 
become the new design point. Otherwise, the design point does not change. The 
search vector multiplier (+/- search magnitude parameter or zero) used with 
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the search vector is stored for later use. This procedure is then repeated 
with the remaining columns of the search matrix. 

A new search vector is created after using all of the search vectors in 
the search matrix. This new vector is created by vectorially adding together 
all of the search vectors times their search vector multipliers. The new 

search vector is a vector sum of previous successful search vectors since 

unsuccessful search vectors have search multipliers of zero. Thus, the new 

search vector represents (stores) the trend of the optimization process. The 
new search vector is scaled such that the magnitude of it’s largest component 

is unity and then is used to replace the first column of the search matrix. 

The procedure is repeated, a new search vector is determined, and then used 
to replace the second column of the search matrix, and so forth until only 
the last column of the search matrix remains untouched. Then an entirely new 
search matrix is created using the random number generator, and the process 
continues. 

If at any time in the iterative process, a new search vector has a 

magnitude of zero (implying all current search directions are not 

beneficial), then a new random search matrix is created immediately. The 
random number generator uses a seed based on the number of seconds f rom 
midnight on the computer’s clock. Each successive run of the optimizer will 

use a different set of search vectors. Currently, the program runs the 
optimizer two times (each time using different sets of random search vectors) 
to help ensure that the global maximum of the coefficient of determination 
has been located in the design space. The number of random search matrices 
generated in each run is equal to twenty times the number of design 
variables. 

The effectiveness of this prediction technique is described in [21. 

6. The Thermal Analysis Program 

A numerical model to predict the thermal behavior of impact damaged MLI 
was developed during this investigation. In this section the theory and 
assumptions associated with the thermal model are discussed. 

The main goal of this project was to develop a microcomputer-based 
design tool to approximately predict the effects of damage to the MLI of 
Space Station Freedom. To be suitable as a design tool 'requires that the 
program be easy to use and that solution times be minimized to rapidly 
provide feedback for design studies. These requirements dictated that the 
numerical model be made as simple as possible while still retaining the 
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capability to provide physically reasonable results. 

The numerical model was based on the assumption of axial symmetry about 
the center of the MLI damage. A finite difference analysis approach was used 
to discretize the system, where an axially symmetric ring of material can be 
approximately modeled as a single node as shown in Fig. 6.1. Higher levels of 
accuracy can be obtained by using more nodes and spacing them closer 
together. Thus, only a single, radial line of calculation points (nodes) was 
required for each layer in the thermal system. The numerical model uses the 
same number of nodes in each layer. The time required to complete a set of 
calculations increases greatly as more nodes are used. 

The computer program is designed to automatically refine the mesh of 
nodes until further refinement produces no change in the results or until the 
user specified maximum number of nodes per layer Is reached. Each refinement 
halves the radial spacing between the nodes. The advantage of this refinement 
process is that a coarse mesh (large node spacing) is used to relatively 
quickly calculate an accurate set of nodal temperatures and heat fluxes which 
are then used as initial values for the refined mesh. Accurate initial values 
for the nodal temperatures and heat fluxes greatly enhance the rate 
convergence to a solution. An accurate solution can usually be obtained 
faster using a series of progressively finer meshes than If a single fine 
mesh is used. Also, the multi-mesh results provide the user with information 
on the sensitivity of the calculated results to the node spacing. 

The pressure wall, the MLI blanket, and the bumper were assumed to 
radially extend out to infinity. The presence of ring frames and stringers is 
not modeled. These would be difficult and computationally expensive to model 
since they would be arbitrarily placed which would destroy the radial 
symmetry. Accurate studies of the effects of the ring frames and stringers 
would require a very detailed special purpose thermal model. Such studies are 
beyond the scope of the design tool under development in this study. However, 
the presence of the ring frame and stringers was accounted for indirectly 
during the thermal model parameter calibration process as discussed in 
13 ]. 

As was discussed, it was assumed that all the MLI consisted of the same 
number of layers, and that no lap joints were present in the MLI. It was 
assumed that the damage in the MLI consisted of a circular hole of the same 
diameter through all of the MLI layers. Deviations from this assumed ideal 
hole geometry are provided by an experimentally determined parameter called 
the "diameter ratio" which will be described. Each layer of the MLI is 
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explicitly modeled with an array of nodes. All aluminized layers are modeled 
in the same fashion. However, the beta cloth layer and the outer kapton layer 
are modeled as a single layer since they are not separated by a dacron 

netting spacer. Thus, for the baseline insulation system of the Space Station 

shown in Fig. 6.2, 22 layers have to be modeled: the pressure wall, 19 
aluminized ML1 layers, the combined kapton - beta cloth layer, and the bumper 
layer. 

The numerical model was designed to model steady state conditions. 
Steady state means that the heat flux into each node In the system must equal 
the heat flux out of that node. Thus, an equation can be written for each 
node to calculate the nodal temperature such that heat influx will equal heat 
outflow. The MLI is in a vacuum so the modes of heat transfer are by 
conduction and radiation. Since radiation heat transfer is a function of 
temperature to the fourth power, the nodal heat flux equilibrium equations 
are nonlinear and must be solved by an iterative process. The thermal 

equilibrium equations are coupled as well - the temperature of a given node 

depends on the temperature of adjacent nodes in the same layer as well the 
nodal temperatures of adjacent layers. This complex pattern of heat flow is 
schematically illustrated in Fig. 6.3. 

As adapted from [10] the following equation can be written to describe 
thermal equilibrium at a node: 


r dr 2 


dr 


q - q ■ q 

In out i 


( 6 . 1 ) 


where: k is the in-plane thermal conductivity of the layer, X is the 
thickness of the layer, r is the radial position of the node, T is 
temperature in the layer, q is the heat flux into the layer at position r 

In 

from adjacent layers, q Is the heat flux out of the layer at position r to 

out 

the adjacent layers, and q ( is the net flux into the i-th node. Eq. 6.1 
basically states that the heat conducted away from a node in the plane of the 
layer must equal the net influx of heat to that node from adjacent layers. 

A standard finite difference approach was used to calculate the 
temperature derivatives at the i-th node: 
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where: A is the radial distance between nodes, is the temperature of the 

i-th node, and T and T are the temperatures of the nodes on each side 
l-i 1+1 

of the i-th node in the layer. The (i-I)-th node is closer to the origin of 

the coordinate system than the (i+l)-th node. 

Note that the (1/r) factor prevents Eq. 6.1 from being used to calculate 
nodal temperatures at the origin of the coordinate system (r = 0 at node 1) 
for the pressure wall layer. The same problem occurs for the special case 
where there is no hole in the ML1, thus the first node of each MLI layer and 

the bumper is at r = 0. This singularity problem was solved in conjunction 

with treating the boundary conditions. For the case of a layer with no hole 
(node one at r * 0) axial symmetry dictates that the in-piane radial heat 
flux through the origin must be zero. This can be ensured by setting (dT/dr) 
and (d^/dr 2 ) equal to zero at node 1. Considering the form of Eqs. 6.2 and 
6.3, this required setting T » T ■ T and so the temperatures at nodes 1 

and 2 were simply set equal to that of node 3. The same approach was also 

used for the MLI layers for the case where there was a hole in the MLI, since 
here there was also no radial flux at node 1 because of the presence of the 
free edge. 

The technique used to treat the boundary conditions at the outer edge of 
the modeled area will now be discussed. The user of the program specifies the 
radius of the area to be modeled and the number of nodes, N, per layer. The 
N-th node would be located at on the outer edge of the modeled area. To 

preserve a type of symmetry in the matrix of governing equations, the 

computer program automatically adds an (N+l)-th node to each layer. It was 

assumed that at the N-th node the perturbing effects of the MLI hole have 

died out. Thus* the radial heat flux would be negligible and the radial 

temperature profile uniform at the N-th node of each layer. This boundary can 

be modeled by setting T and T equal to T . This same boundary condition 

would apply if the boundary of the area modeled was aligned with the outer 
edge of the pressure wall plate, the bumper plate, and the MLI blanket. Thus, 
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the inner and outer boundary conditions for every layer were treated in an 
identical fashion. 

Substituting Eqs. 6.2 and 6.3 into Eq. 6.1 produces the following 
equation which describes thermal equilibrium at the 1-th node of a layer: 
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Eq. 6.4 can be written in a more compact form as: 


(6.4) 
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C T - q 
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(6.5) 


where the C (j = 1 to 3) are thermal equilibrium influence coefficients for 
the i-th node of a layer which can be evaluated from Eq. 6.4. Eq. 6.5 can be 
expanded in matrix fashion to represent an entire layer of nodal 
temperatures: 
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Eq. 6.6 consists of a tridiagonal system of equations which is very efficient 
to solve numerically. Note that T and T are not explicitly solved for, 
rather they are set equal to T and T , respectively, of the previous 

3 N-l 

iteration as was mentioned in the boundary condition discussion. An iterative 
procedure is required to solve Eq. 6.6 because the qj values are complicated 
nonlinear functions of the nodal temperatures. 

The solution procedure consisted of solving Eq. 6.6 for the nodal 
temperatures of the first layer (pressure wall) and then proceeding to the 
next layer, solving for the nodal temperatures and so forth, until finally 
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solving for the nodal temperatures of the final layer (bumper). This 
procedure is repeated until the nodal temperatures converge with respect to a 
user defined tolerance. 

The calculation of the q^ values will now be discussed. The formulas 

used to calculate the values varied from layer to layer. Accordingly, the 

method of q calculation will be discussed on this basis, 
n 

PRESSURE WALL NODAL NET HEAT INFLUX CALCULATIONS 

The pressure wall will interact thermally with the atmosphere of the 
inside of the spacecraft In the form of forced convection heat transfer. The 
equation describing this process is 1 10]: 


q - h(T - T) (6.7) 

C 00 | 

where q^ is the convected heat flux to the i-th node of the pressure wall, h 

is the convective heat transfer coefficient, T is the free stream air 

00 

temperature of the spacecraft module, and T is the temperature of the i-th 
pressure wall node. Thus, if the pressure wall is cooler than the module air 
temperature, then heat from the spacecraft module air will flow into the 
pressure wall and vice versa. 

The convective heat transfer coefficient can be estimated from the 
following equation which was adapted from information presented In [10]: 


h « 4 



1/2 


( 6 . 8 ) 


where u is the velocity of the module air next to the pressure wall (m/s) 
00 

and L is the distance (m) that the air travels along the pressure wall before 
meeting an obstruction such as a ring frame. 

The pressure wall radiates heat towards the MLI blanket and into the air 
of the spacecraft module. This heat flux, q^, is described by the following 
equation [10]: 


q r = co-tJ (6.91 

where e is the emissivity of the radiating surface, <r is the Stef an-Boltzmann 
constant [5.6697E-8 W/(m 2 K 4 )], and T is the temperature of the radiating 
node. Emissivity values can vary between 0 and 1 depending on the material 
the surface is made from and the condition of the surface (polished or 
tarnished and so forth). Emissivity values can vary as a function of time and 
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temperature. In this investigation all emissivities were assumed to be 
constant. 

The pressure wall was exposed to heat flux radiated down from the 
adjacent MLI layer, and from the bumper and space environment if a MLI hole 
is present. Not all of the radiation impinging on the pressure wall was 
absorbed. The fraction absorbed is called the absorptivity. To simplify 
calculations, the absorptivity is commonly assumed to equal the emissivity 
[10]. Radiated energy that is not absorbed is reflected. To preserve 

conservation of energy, the computer model keeps track of the magnitudes of 
emitted, absorbed, and reflected radiation. Actually, a portion of the 
radiation striking the pressure wall from the adjacent MLI layer and from the 
bumper and space environment through the MLI hole is reflected radiation from 
these layers. 

The simplest case of no MLI hole will be considered first. Here, all the 
nodal temperatures of the MLI layer next to the pressure wall will be 

identical after equilibrium is attained. Thus, the thermal radiation emitted 
and reflected will be the same for each node in the MLI layer. Also, no 
thermal energy from the bumper or space environment will strike the pressure 
wall. Accordingly, for this simple case, the computer program uses the 

thermal radiation (both emitted and reflected) from i-th node of the MLI next 
to the pressure wall when calculating the heat influx to the i-th node of the 
pressure wall. 

The more general case with a hole In the MLI i9 considerably more 

complicated. Here, the thermal radiation coming from each node of the MLI 
layer next to the pressure wall will vary. Also, the thermal radiation from 
the bumper and space environment will pass through the hole In the MLI and 
strike the pressure wall plate. The concept of view factors [10] was used to 
treat this problem. 

View factors give the fraction of the thermal radiation given off from a 
surface that will strike another surface of known geometry and position. 
Consider Fig. 6.4 where thermal energy Is radiating from circular area A^ to 
circular area In Fig. 6.4, the plane of area A^ is parallel to the plane 
of area A^. The view factor associated with this geometry, F^ ^ is given 
by [10]: 
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1 + G 2 + H 2 - /(I + G 2 + H 2 ) 2 


( 6 . 10 ) 
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AI-A2 


2 2 
4G H 


2G 


where G = b/a and H * c/a (see Fig. 6.4). Note, for example, that F 

AI-A2 

approaches unity as c (and thus G) approaches infinity as one would expect 
because for this case area A ^ would be radiating into an infinite plane and 
thus all radiation would be captured. 

As is illustrated in Fig. 6.1, the numerical model developed during the 

course of this investigation is based on the assumption of axial symmetry. 

Thus, a view factor, F , for radiating from ring area to ring area is 

r 

required here. F can be obtained by repeatedly applying Eq. 6.10 (see Fig. 

r 

6.5): 
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A (F - F ) - A (F - F ) 

11 A 1 1- A21 A 1 1 - A 2 2 12 A12-A21 A12-A22 


1 2 


( 6 . 11 ) 


Thus, F specifies the fraction of the energy that is radiated by AA (= A - 

r 1 11 

A ) that will strike AA (=* A - A ). 

12 2 21 22 

The total heat influx to the ring corresponding to the i-th node of the 

pressure wall from the adjacent MLI layer was calculated by a summation 

formed from repeatedly using Eq. 6.11 for each node (and thus corresponding 

ring) of the MLI layer. This is shown schematically in Fig. 6.6. The outer 

boundary of the ring corresponding to the N-th node of the MLI layer was 

extended out a large distance beyond the user specified radius of modeled 

area. This was done to be compatible with the assumption that the layers 

extend out to infinity in all directions. 

The heat Flux from the bumper to the pressure wall through the hole in 

the MLI was treated in two steps. First the ring approach of Eq. 6.11 was 

used to calculate to total heat flux to the MLI hole from each node on the 

bumper. For this calculation, A (Fig. 6.5) was set to zero and A 

22 21 

corresponded to the area of the MLI hole. Also included in this calculation 
is the thermal radiation from the space environment (space thermal radiation) 
which would pass through the bumper hole. Then, the thermal energy impinging 
on the MLI hole was allocated to the pressure wail node under consideration 
by using Eq. 6.11 again. For this calculation, A j2 was set to zero and A^ 
was set equal to the area of the MLI hole. This process Is Illustrated in Fig 
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6,7. 


The sun continuously emits thermal radiation, a small fraction of which 

strikes the earth. This incident solar radiation flux from the sun has an 

average magnitude (called solar constant GJ of approximately 1353 W/m [10]. 

The thermal analysis program uses a parameter (in file THERMAL.PAR) called 

space thermal radiation to account for thermal radiation striking the bumper 

from the space environment. Thus, if the region of the spacecraft wall being 

modeled directly faces the sun, then the space thermal radiation parameter 

should be set to 1353 W/m 2 . If the modeled region faces deep space then this 

2 

parameter should be set to 0 W/m . A cylindrical module with one side facing 
the sun and the other facing deep space would experience an average heat flux 
of 431 W/m 2 . 

NODAL NET HEAT INFLUX CALCULATIONS FOR MLI LAYER NEXT TO PRESSURE WALL 
The MLI layer next to the pressure wall (first MLI layer) can radiate 

energy to both the pressure wall and the next MLI layer. Thus, the q for 

r 

this layer will be twice that given by Eq. 6.9. 

The pressure wall can subject the nearest MLI layer to both emitted and 
reflected thermal radiation. This was treated in exactly the same way that 
MLI heat flux impinging on the pressure wall was treated (reverse of Fig. 

6.6), which has been discussed. Note that the MLI layer next to the pressure 
wall Is blocked from receiving radiation directly from the bumper or the 
space environment. 

The MLI layer nearest the pressure wall will also be subjected to 

emitted and reflected radiation from the next MLI layer (second MLI layer). 

Since the MLI layers are so close to each other a view factor approach of Eq. 

6.11 was not used here. The thermal radiation flux from the second MLI layer 

striking the i-th node of the- first MLI layer was assumed to equal the 

thermal radiation flux from the I-th node of the second MLI layer. 

Direct conduction between the first and second MLI layers was inhibited 

by the presence of a layer of dacron netting. The heat flux to the first MLI 

layer from the second MLI layer through the dacron netting, q , was assumed 

N 

to be of the following form: 


h (T - T ) 

N 1,2 1,1 


( 6 . 12 ) 


where h is the effective netting heat transfer coefficient, and T and 

n 6 1,1 

T are the temperatures of the i-th node of the first and second MLI 

1,2 
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layers, respectively. A value for h equal to 1.0687 W/m Z K was determined for 

N 

the Space Station MLI by fitting the computer model to experimental data as 
discussed in [3]. It was assumed that the netting heat transfer coefficient 

was the same for all netting layers. 

NODAL NET HEAT INFLUX CALCULATIONS FOR A TYPICAL ALUMINIZED MLI LAYER 
Here the net heat influx, q , to the nodes of the MLI layers between the 
the first (next to pressure wall) and last (next to bumper) MLI layers are 
considered. The q ( values for the nodes of these layers are calculated in a 
similar fashion to what was done for the first MLI layer, except here there 
are two MLI layers radiating into the MLI layer under consideration. No view 
factor calculations are required here, since the layers are assumed to be 
close together. Also, there are two layers of dacron netting next to each MLI 

layer, and thus Eq. 6.12 will have to be applied twice - once for the layer 

above and once for the layer below. 

NODAL NET HEAT INFLUX CALCULATIONS FOR LAYER NEXT TO BUMPER 

As was noted previously, the last aluminized MLI layer (closest to 
bumper) and the beta cloth layer are not separated by a layer of dacron 

netting (see Fig. 6.2). Accordingly, these layers were analyzed as a single 
layer with the inside surface having the emisslvlty of an aluminized layer 
and the outside surface having the emissivity of the beta cloth layer. The 
thermal conductivity of the layer was assumed to equal the weighted average 
(on the basis of thickness) of the two layers. For q^ calculation purposes, 
this combined layer was treated In exactly the same manner as the first MLI 
layer except that here the takes the place of the pressure wall, 

NODAL NET HEAT INFLUX CALCULATIONS FOR THE BUMPER LAYER 

The net heat influx to the bumper layer was calculated in a very similar 
manner to that of the pressure wall. The bumper will be subjected to heat 
influx from the pressure wall through the MLI hole Just as the pressure wall 
was from the bumper. However, unlike the pressure wall, there is no 
convective heat transfer to the bumper from the spacecraft module air. 
Instead, the bumper interacts with the thermal radiation from space. 

This concludes the discussion of q ( calculation for the various layers 
of the thermal system. 

Nodal temperatures were calculated layer by layer starting with the 
pressure wall, and finishing with the bumper. A set of calculations covering 
all layers once is considered one global iteration. The user sets the number 
of global iterations in the thermal parameters file THERMAL.PAR. 

It typically takes many global iterations until the nodal temperatures 
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converge. After convergence is reached, the program refines the mesh and then 
starts calculations for the new mesh. If the finest mesh is being used when 
global convergence occurs, then the program stops. If the maximum allowable 
number of global iterations (a user input parameter) is used before 

convergence is obtained, then the program refines the mesh and begins 
calculations again. If the finest mesh is being used and convergence is not 
obtained before the maximum allowable number of global iterations has been 
exhausted, then the program issues a warning and stops. 

Ten global iterations are conducted between each check for convergence. 
Convergence Is assessed by calculating magnitude of the change that occurred 
in the temperatures of the Inside and outside edge nodes of the pressure wall 
and bumper layers during the ten global iterations. The change in temperature 
is divided by the magnitudes of the temperatures to produce a nondimensional 
relative temperature change. The relative temperature change is compared with 
a user input convergence factor stored in the thermal parameters file. If the 
relative temperature change is less than the convergence factor, then the 
calculations were considered to have converged. 

As has been discussed, the computer program has been designed to 
automatically refine the mesh by halving the distance between the nodes. The 
idea is to have coarse meshes provide accurate initial values for 

successively finer meshes. This serves two purposes: the rate of convergence 
is enhanced and Information on the sensitivity of the calculated results to 
the mesh density is provided. Ideally, the mesh should be refined until there 
is an acceptably small change in the calculated results. 

A factor called "diameter ratio” was developed In 13] to account for the 
fact that the apparent MLI hole diameter with respect to thermal behavior 

tends to be different from the measured MLI hole diameter. Thus, the diameter 
ratio accounts for damage effects such as the charring and crinkling of the 
MLI beyond the edge of the MLI hole. It also indirectly accounts for the fact 
that the MLI blanket has some thickness. The diameter ratio is defined as the 
thermally apparent diameter ratio divided by the visually measured diameter 

ratio. Thus, for thermal calculation purposes, the visually measured MLI hole 
diameter should be multiplied by the diameter ratio parameter. 

The empirical function for diameter ratio that was derived during the 

course of the investigation reported in [3] is: 


31 



diameter ratio = D 0.1978 V 3 ‘ 466 T 5 ' 356 D _7135 (cos 0) 1 ’ 694 + 2.575 

P >P 

(6.13) 

where D is the diameter of the projectile, V is the velocity of the 

p 

projectile, T is the thickness of the bumper, and is the impact angle (see 

b 

Fig. 1.1). Of course, Eq. 6.13 should only be used to predict diameter ratios 
for impact conditions similar to those of the investigation reported in [3], 
Otherwise, a diameter ratio of unity can be used as an approximation. The 
form of Eq. 6.13 was derived from the nondimensional prediction equations 
discussed in Section 5. 


32 



7. The Condensate Prediction Program 

7.1 INTRODUCTION 

The condensation process of water vapor from moist air over a circular 
surface is studied here. The main objective of the study is to determine the 
condensate height for a given temperature distribution on the surface. 

Typically, condensation problems have been dealt with using boundary layer 
techniques. Two sets of conservation equations are solved: one for the 

condensate layer and one for the vapor layer, with appropriate interface 

conditions. But the boundary layer theory breaks down near the center of the 
circular region [15]. Hence the full Navier-Stokes (momentum conservation) 

equations are to be considered. 

In the following case, the three basic conservation equations of mass, 
momentum, and energy are solved with changes in thermophysical properties 

being accounted for with changes in temperature. The height of the condensate 
is determined from the final temperature distribution over the region. The 
input parameters required are the radial positions of the nodes along the 

pressure wall and the corresponding temperatures, the radius of the surface 
on the pressure wall, the ambient and dew point temperatures of moist air, 

and the velocity of air over the surface. 

7.2 LITERATURE REVIEW 

Analyses of laminar film condensation problems have generally been done 
for flow of a vapor or vapor-gas mixture along horizontal or vertical 
surfaces [18-251. All of the studies dealing with vertical surfaces involve 

gravity as a body force. The case of film condensation on a horizontal fiat 

plate in the absence of gravity has been studied [25]. But in this study, a 

boundary layer formulation has been used and as explained earlier, this 

solution breaks down in the vicinity of r = 0. It should be noted that in 

most of the above cases an uniform wall temperature was assumed. 

7.3 PHYSICAL MODEL AND COORDINATES 

The model for the problem along with the coordinate system is shown in 
Figure 7.1. The flow of moist air Is directed radially away from the center 
with a velocity v in j et - F or an axisymmetric flow situation, the centerline 

r = 0 , which forms the left boundary of the domain, is an axis of symmetry. 

The wall along which the temperature Is prescribed as a function of the 
radial position, T * T(r), forms the front boundary. Since the edge of the 
circular region is Insulated, the temperature equals that of the ambient 

condition for r t rmax. The right boundary (r » rmax) is an open outflow 
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boundary. For the present case, it Is assumed that ambient conditions are 
achieved at a distance of one radius from the surface in the axial direction, 
i.e. zmax ■ rmax. The condensate layer will extend in the radial direction 

till the point at which the wail temperature exceeds the dew point 
temperature, while the height in the axial direction will vary with r 

depending on the value of z at which the vapor reaches its dew point 

temperature. 

7.4 GOVERNING EQUATIONS 

The general form of the conservation or transport equation in 

cylindrical coordinates (for an axlsymmetric problem) for any variable 9 is 
given by U6) 


3(p*) 1 d(pru4> ) d(pv<t>) 1 3(rr fj) 3(r |S) 

St 7 87 dz “ 7 S7 * 3r * 55 * dz + s * (7 - 


a*. 


9 


1 ) 


The first term on the left is the unsteady term which accounts for changes In 

9 with respect to time. The second and third terms on the left represent 

changes In 9 due to convection. 9 Is called the convected quantity while u 

and v are the convectlng velocities. The first two terms on the right 

represent changes in 9 due to diffusion, r, is called the diffusion or 

9 

exchange coefficient. The last term on the right is called the source term. 

Effects not represented in the other terms, for example the pressure gradient 

term, are included in S,. 

9 

For 9 = 1. eqn.(7.1) becomes the mass conservation or continuity 

equation. If 9 is taken as u or v, we get the momentum conservation or 

Navfer-Stokes equations. If 9 is the enthalpy h, we get the energy 

conservation equation. Since h ■ C T, where C is the specific heat at 

P P' 

constant pressure and T is the temperature, the energy conservation equation 
can be written in terms of T. The conservation equations mainly differ in the 

form of the diffusion or exchange coefficient and source term S, as shown in 

9 

Table 7.1. Here ji is the absolute viscosity, k is the thermal conductivity, 
and j s the substantial derivative and Is written as 

D3 8 8 

Dt at + u a? + v ai 

For an incompressible flow situation, p can be considered to be a 
constant. Hence the continuity equation Is 


l a(ru) 8v 

r dr dz 


(7.2) 


34 



The momentum conservation equations are written as 
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where i> = ^ Is the kinematic viscosity. 

For the energy conservation equation, we write 
Considering C ^ to be constant we have 


the enthalpy h= 


C T. 

P 


ar i a(ruD a(vT) 
at + r 3r + 3z 
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1 Dp 
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(7.5) 


where a is the thermal dtffusivity which equals 



7.5 COMPUTATIONAL PROCEDURE 
7.5.1 Description Of The Grid 

The computer code for solving the governing equations Is based on the 
marker and cell (MAC) method. This method uses a nonuniform mesh system. The 
scalar quantities like pressure, temperature, and density are placed at the 
center of the cells while the velocity components are placed normal to the 
cell walls as shown in Figure 7.2. This kind of placement of variables Is 
called a staggered mesh system. Hence, there are three different control 
volumes or cells for u, v, and T which are used in solving the momentum and 

energy conservation equations, Figure 7.3. The boundaries for the cells are 

denoted as the left, right, front, and aft faces. 

A fictitious layer of cells is added to all four sides of the 
computational domain. Hence, at any of the boundaries the normal velocity 
components lie directly on the boundary while the tangential velocity 
components and any scalar quantities are displaced by half a cell width 

within the flow domain. The additional layer of cells makes it easy to apply 
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various types of boundary conditions as shown in the next section. 

The widths of the T-cell are Arj and A Zj in the radial and axial 
directions respectively. Since the velocities are displaced half a cell 
width, the u and v cells used in the computation are displaced by the same 
distance in the radial and axial directions respectively. The widths of the 
u-cell and v-cell in the r and z directions respectively are defined as 


Au, 


Ar^/2 + Ar l+1 / 2 


Av 


A z j/Z + LZj^/2 


The different dimensions of the cells are shown in Figure 7.3. The distances 
to the center of the cells are rc^ and zCj while the distances to the u and v 
locations are ru^ and zv^ in the radial and axial directions respectively. 

7.5.2 Boundary Conditions 

It is assumed that ambient conditions are reached at a distance equal 
one radius from the origin in the axial direction. Accordingly, the dimension 
of the computational domain in the z-direction is taken to be equal to the 
radius. 

The moist air enters the domain through the aft boundary, i.e, flow Is 
radially away from the center. The front boundary is a rigid no-slip surface 
and the velocity components are zero on this boundary. The temperature is set 
according to the given distribution. The left boundary of the domain is the 
centerline (r*0) and hence a line of symmetry. Thus, the normal velocity 
component, u, is zero on this boundary while the tangential component, v, and 
the temperature have zero normal gradient across the boundary. The 
right boundary of the domain is the outflow boundary and continuative 
conditions are applied here. The gradients for both the velocity components 
are set equal to zero. The temperature along the right boundary was set equal 
to the ambient temperature. The following boundary condition was imposed at 
the aft (inlet) boundary. The axial velocity v was set equal to the inlet 
velocity v j n j et and the radial velocity u was set equal to zero. 

The boundary conditions can be summarized as follows: 

Aft Boundary (inlet) 


L t Jmax 

i tjmax 


i t Jmax-\ 

~ U l t Jmax-l 


inlet 
(u « 0) 


( T + T )/2 ■ T 

u l,Jmax l.Jmax-l 1 ambient 
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or 
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Front Boundary (rigid no-slip) 
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where T(r) Is the given radial temperature distribution. 


Left Boundary (symmetry) 
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-u 


2 J 


(u - 0) 


V 1J “ V 2J 

T 1 ,J " T 2J 


Right Boundary (outflow) 

u imax,J * U lmax-l,J 

V my 

lmax,J lmax-\.,J 

T ■ 2(T ) — T 

1 lmax,J ambient' lmax-l,J 

The outflow boundary conditions are imposed only after each time step and not 
after each step of the pressure iteration procedure which is explained later 
on in the text. This is because the normal velocity at the outflow boundary 
may vary with changes in pressure. 

7.5.3 Solution Method 

The steady state problem of condensation in this case is solved using an 
unsteady time-dependent technique [16]. The solution proceeds by marching 
forward in time, with the time-marching procedure being continued till there 
is negligible difference in the values of the variables between two 
consecutive time steps, i.e. a steady state has been reached. 

Since it is an incompressible flow problem, the density is not 
considered to be a variable but is updated after each time step since it Is 
temperature dependent. The other transport properties like viscosity and 
thermal dlffusivity are also updated after each time step. 

The governing equations are solved by means of an explicit finite 
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difference method [16]. The subscripts i and J denote the radial and axial 
directions respectively while the superscript n stands for time level t. No 
superscript is used for time level t+At. The detailed derivation of the 
finite difference approximations to the governing equations is given in the 

next section. The convection terms are discretized using the first upwind 

differencing scheme while the diffusion terms are evaluated using central 

differencing. Forward differencing is used for the time derivatives. 

At each time step the velocities u and v are calculated explicitly from 
the momentum conservation equations and the temperature T from the energy 
conservation equation. The calculated velocities are used as initial values 
for the next time step. These velocities do not in general satisfy the 

continuity equation. .The reason for this is that the pressure field is not 

known a priori (an initial guess for the pressure field has to be made at the 
beginning of the solution) and since velocity is affected by pressure 

changes, the cell pressures have to be adjusted such that the velocities u 
and v satisfy the mass conservation equation. This procedure is explained 

later on In the text. 

7.5.4 Finite Difference Approximations To Governing Equations 

Following are the finite difference schemes used to evaluate the 

different terms of the governing equations: 

(I) Time Derivatives - Forward Differencing 

(II) Convection Terms - First Upwind Differencing Scheme 

(iii) Diffusion Terms - Central Differencing 

(iv) Source Terms - Central Differencing 

Though central differencing could be used for the convection terms too, 
it has been found that [26] representing the convected quantity ^ by central 
differencing leads to instabilities in the solution. Let us take the simplest 
example of an uniform one-dimensional mesh with x as the coordinate and a 
uniform velocity u through the mesh. 

*i-i m 4 i *r 8 i 5 

o i O , O > X 

(-1 L! i !r i+i 

Mesh for one dimensional problem 

„ .j . S(urf>) 

Consider the term ^ . Using central differencing, we have 

*R " ( *i + *i+l )/2 : * l *i + 
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a(_u<t>) a U R* *R ~ U L* *L = ( *i + *t + l )/2 " ( *t + Vl )/2 

9x ~ Ax f U AXj 

^i + 1 " ^l-l 

= u - (since u R = u L = u) 

This Implies that changes inside the l** 1 cell are not affected by the value 
of but only by those of and It Is clear that whatever the 

values of <f> and may be, the value of 4 > could still lie outside this 

range. For example, let 4, <f> 8, and 4 5 as shown in the mesh 

diagram above. Assuming that no sources or sinks are present, the <f> field has 
to be monotonically increasing or decreasing and no sudden jumps are allowed. 
Hence the above field is not physically plausible. 

The easiest alternative to this is the upwind differencing scheme [26] 
in which the value of 4> at any cell face is chosen to be the upstream value 
according to the sign of u at that cell face, i.e. 

*R = +1 ‘ f U R > ° . *L = *£- 1 If U L > ° 

I 

= f , if u. < 0 = <p. if u. < 0 

l+l R i L 


All other quantities whose values are not directly available at the cell 
faces are obtained by linear interpolation from the adjacent cell values. The 
thermal conductivity k is represented in a different manner than the other 
properties as explained below. Again, let us consider the uniform one 
dimensional mesh. The heat flux at the right face of the cell is given by 


q R 




(+1 


Now consider the cell between the t** 1 and 1+ 1 th grid points as a composite 
slab. Then the heat flux, from basic principles, is 


i 


i + 1 


Ar { /2 Ar l+1 /2 

k, + )T 


l 


l+l 


Comparing the two expressions for q R , we get 
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Au 


Ar (+J /2 Ar ( /2 

— : + — T 


R 


( + 1 


or 


_ k l * *1 + 1 
*R _ k.* A r 


1 

* A u 


l 


l 


i./ 2 * *i.i* ir i /2 


We use this form since it gives the exact representation of the heat flux 
across the cell face. 

The r-momentum equation is discretized with respect to the u-cell. 


n 

(i) au „ U j,J ~ u i,J 
at ~ At 


(ii) . . 1 fl(ruu) 1 rc i + l* U R?u* U R " rc t* U L?U* U L 

{a) 7 d? “ FUj AU] 


“ R ?u - (u <+l 


+ 1 U M =(U tV U KJ ,/2 


~ n .r n ^ r\ 

U R ‘ “l.J lf U R.u > 0 

= u ", , lf u " < 0 

i+l,J R.u 


~ n ir n . . 

U L ■ V l.J lf U L. U " 0 

= u" . if u " <0 

l.J L.u 


. . v n * u, - v n * u_ 

, L , a(vu) „ A,u A F,u F 

,b) 55 =; 


v, n ,* Ar,.,/2 + v i+i f j* Ar j/2 


n _ _ 1+T 

/ A, u 


Au 


l 


v" . * Ar /2 + v ” . .• Ar /2 

n 1 , j-1 t + 1 t + l.J-1 i 


F,u 


Au 


l 


n .. n _ 
u = u. . if v. > 0 
A l,J A,u 

= u, . , if v. < 0 
i A , u 


~ n n 

U F ’ lf V F,U > 0 

= u n . if V_ n < 0 

l,J F , u 
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_ du * ou | 

au ' ■ ri, ^lp - rv a?\ L 


du i 


(m) , y 1 dirv-^) 1 3r R 0r 

(a) — — or « t 

r 3r ru, Au 


t 


i 


n n n n 

u. f . - u. . u. . - u. - . 

rc « v n ( . IhJ. LlA) - rc * v n ( t—hJ.) 

rc l + 1 V R,iT Ar,., ' 1 L,u l ir, 1 


l+l 


i 


ru l* Au t 


n n n n 

v R,u ~ V l+lJ : V L,u = V i,J 


_ du | mm 0u * 

fh) a(^) „ azU dz I f 
(b) ai 5 IT. 


n n n n 

n 1 ' “l.J> .. n , U l.J " U l,J - 1 

AV J " F ‘ U &V J -1 


f » / 


) 


Az 


Values of p n and p_ n are obtained by linearly interpolating from the four 
A f Li r , u 

neighboring v values. 


A,a 4Au 


J 


v A r Az 
i,j+r l+l 


vJ 1 * — — fy. n , ,Ar.Az . + p. n , . Ar Az . + v n .Ar Az 

F,u 4AUjAv^ j[ l+l (v r i J-l l+l.j-1 i J l,J l+l J-l 

■J 


l ’i,j-i lr i.i iz 
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(iv) (a) 


The 


n n 

l iE » - L il±Ll ~ p _Ld 

p dr p A u . 

^ u L 


u 


(p ( . Ar {+] /2 + p {+1 * Ar t /2) 
Au . 


„ 3v i 9v . 

|j.^) . V|a_ V>f 

9z Az . 


n (. 1 LiJ.) - v " ( 

A.u 1 Au ( F,u 

Az , 


n n 

n , tfl.J-1 ~ Vh, 


Au 


i 





ru 


2 

( 


v 


n 

u 


- iv 


n 

t+W 


+ 




)/2 


z-momentum equation Is discretized with respect to the v-cell. 


n 

(i) Vj 

at At 


1 a(ruv) „ 1 ru t* u R?v* V R " ru t-l* U L?v’ V L 


(ii) (a) 

r dr rc 


i 


hr 


n 


u" . Az , ,/2 + u, n , .• Az ./2 


U D - Jbi - 

R , v Av 


t.J+1 J 


n _ u l-l,J m Az J+l /2 * u i-l, J+l* Az J /2 


L,v 


Av 


* 
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" Vv > 0 

= v. , , if u n <0 
i+lj R,v 


v ■ v n if u n >0 
L l-l ,J L, v 

- v n if u. n <0 
i.J L, v 


(b ) £i v ' v ') * V A?v* V A ~ V F?v* V F 

dz Av , 


V A,v " (v i,J+l + V t,J )/Z : V F?v " (v t?j + V 1 ?J- 1 )/2 


v = v" if v n > 0 
A l,J A, v 

n n _ _ 

= v. . . if v. < 0 
l,J+ 1 A, v 


v = v n if v n > 0 
F t.j -1 F , v 

« v n ir v n < o 

t.J F , v 


, , 1 d{rv~) 1 

(,,1) (a) r Jr dr 


dV • uv | 

rv s? Ir - '"irk 

A r. 


3 v i 


t 


n n n n 

n . v t+l, J * V l.J, n , v l,J " v t-l,J. 

ru.* ( 4 — ) - ru. ,• v , ( — — ) 

t R,v Au. t -1 L,v Au. 


I 


t-1 


rCj* Ar^ 


The values of v n and v. n are obtained by linearly interpolating from the 

K , V L| V 

four neighboring v values. 


n 


v 


R » v 



l ’l"l,J 4 r i 42 Jtl] 

44u ( | Av j ["l , >l 4r l-l Az J * 1 'l-l,J+l 4r l 4 *J * l 'l./ r i-i 1Z J*l * 
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(b) 


dv i 


3v, 


9z * A 


Av 


n 


„ " £hj*i ~ 

A,v l Az J+J 


, - p " ( iLilllLtl) 

F , v Az , 


Av 


n n n n 

A.v * P (,>1 ; ^F, v = 


(iv) (a) 


n n 

- 1 *£ * . 1 P t.>i " P t.J 


p 3z 


Av 


p (?/ Az / + 1 /Z + p l?J + 1* Az / 2 

Av , 


o au | ou « 

1 d(n,i“) 1 ri, ailR - rv al I 

?3F az “7^ 3? 


ri V Vv 1 


n n 

n , U t,J+l ~ U t,J 


Av 


n n 

, n , u t-l,J+l “ u l-i , J 

) - ru. .• v. ( j- ^ 

t-1 L, v Av . 


rCj* A r { 


The energy conservation equation Is discretized with respect to the 7-cell. 


7 - 7 

(i) £I«_Li AA 

at At 


I fUr-,,Ty 1 rU l* U D X* 7 " rU t I* U I T* T , 

(iD ( a ) I £1 ru7) » _i 1 — HiZ — £ L£ Li Z L 

r dr rc^ Ar 


n n n n 

U R,T * u iJ ; U L,T " VlJ 
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f R " T IJ lf “R ?T > ° 


f L S T "lj lf U L?7 > ° 


C.j lf Vr <0 


T C?J ,f U L?T < ° 


, a(vT) V A ,T* " V F , T* 

(b) 


n n 

V A,T = V l t J 


n n 

V * v 

F,T V l t J - 1 


? A - 7 "j ,f V A?T > 0 


f =T n if v n > 0 
F F,T 


T i?Ja. If -*?T ‘ 0 


T t n .j " v f"t < o 


_ dT , ar, 

..... , , 1 3(r«|I) _ 1 ret 3r I R ' ra 3r II 

(in) (a) rTr Sr “ rc^ SF^ 


n m m n ^ ^ 

r V tt R?T ( ltl, tu | J '' ) - rU t-l* a L?T ( t,J Au ~ 


t 


r V 4r i 


i-i 


a 


n n . 
a. .* a. , ,• Au 4 

n = t,J t + l,J i 

R * T ri . n « 

«!./ ar Ul / 2 *“i-l,j* 4r / 2 


a, 


*t? i* au .-i 


L ' T «,?/ 4 W 2 * - n 


1-1./ ir l /2 
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dT, dT, 

, M „ a 3z' A 3z If 

(bl Si dz Si 


'j 


yi n yi ^ n pi i n 

a n ( _Ll*Li1_I l±l) _ a " (J^2___Liz!) 

a.t' Av^ 1 f,t v Av j-i 

Az , 


a 


<V “iV Av j 


A,T n . n . 

*i,/ 4i! j.i /2 * 42 / 2 


n n 

a. .* a. . Av. _ 

, n = l_J t,J-l J - 1 

F,T n . n A ' _ 

*«./ Lz j-i /z + Az / 2 


(lv) (a) 


at 


At 


, dp P R,T " P L,T 
(b) U 5r U T AT; 


u t - °- 5 * (u i,j + u i-i.y 


R,T 


p t,/ Ar t+i /2 * P t+1,/ A V 2 

Au, 


L,T 


p t,/ Ar t- 1 /2 + w Ar t /2 


Au 


£ - 1 


. . 3p P A,T ” P F,T 
M V 3 i * V T Az_, 


V?. - 0,5.(v lJ * v u _,) 
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A J 


p i./ Lz W 2 l * p u*r a v 2 


Av 


P, ,■ + Pi.., -I* to / 2 

F - T 4 'j-l 


The mass conservation equation is discretized with respect to the T-cell. 


1 8(rul .. 1 ru l R ~ rU I L 

r dr rc l Ar ( 


ru * R = r V U IJ ! ru l L " r Vl* VlJ 


(b) 


V - v 

flv w t.j c, J-1 

3z Lz j 


7.5.5 Pressure Iteration Procedure [16] 

The discretized form of the mass conservation equation is given by 


1 (r vi,.< - ru t-iVi.. ) ) , Vj - v i„/-i 


rc 


t 


Ar 


l 


Az 


Since the velocities computed at each time step do not in general satisfy the 
above equation, let the left hand side be denoted as D^(not equal to zero). 
At every time level, the velocities are adjusted by adjusting cell pressures 
until is sufficiently small for all i and J. 

The pressure adjustment is obtained using the momentum equations. 

The discretized r-momentum equation is 


n 


\j 


* i.j 


At P i + 1 ,J P t,J 

P Au, 


+ other terms 


If p. . is Increased by Ap without changing anything else, the adjusted 
velocity is given by 
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n at p l*l.J - p l.J M “ P IJ .. . 

u = u ^-r 41 + t — + other terms 

t.J t.J P Au ( p Au { 


or u, = u, , + 


M Ap t(J 


l,J IJ p A u. 


Similarly u. . , = u 


At Ap 


i.J 


l-l.J “i-l.J p Au 


l-l 


At Ap. 


v, , = v, , + 




t,J l.J P Av 


AtAp^ 


V U - 1 " V iJ-l P Av 


J~ 1 


Now let the adjusted velocities satisfy the continuity equation. Then by 
substituting the expressions for the adjusted velocities in the continuity 
equation, we get 


At Ap 

ruju, , + — — - ru. . ( u 


At Ap 


j t t.J PAu { 


LJ* 


l-l l-l.J pAa 


l-l 


rc 


Ar . 


I l 

At Ap At Ap. . 

+ * »«/ ) _ / .. _ l, J 


(v + - (v 

1 l,J pAv . 1 t.J-1 PAv^_ j 


Az 


= O 


or _!_ (ru i u »../ ~ r “i-i u i-i ..) 1 „ v t. ./ ~ v i..m 

rc i LZj 


. M Ap l,J r 1 r ,u t . fU l-h 1 r 1 1 il _ „ 

+ p IrCjAr-j l-Au ( + Au^jJ + A Zj Iav^ + Av^ ^ J J 


ru t ru 


or 


iP !.J 1 

1 

r ; u ‘ + 1 r * . 1 ii 

P 1 

K ir ! 

LAu Ku J A z Lav Av J 

i i-\ J J J - 1 J 


Denoting the term in braces as 0 


U 


we have 
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(7.6) 



p n 

ITT 


LL 

ij 


It can be seen that the velocity components in each cell are affected by the 
pressure adjustments in all the neighboring cells. So the for all the 

cells need not necessarily be close to zero after just a single step of the 

pressure iteration procedure. Hence the cell pressures and velocities are 

adjusted over and over again in an iterative process, the most recent values 

of velocities being used to calculate D^. The process Is repeated until the 

for all the cells is less than a specified small number e. The right hand 
side of eqn.(7.6) is multiplied by a relaxation factor w to enhance the 

convergence of the pressure Iteration procedure. The optimum value of the 

relaxation factor is approximately w Q p t K 1*8* 

The temperature is computed from the energy conservation equation using 
the velocities satisfying the continuity equation which are obtained after 

the pressure iteration procedure. 

7.6 STABILITY CRITERION 

The time step At used in the computation cannot be greater than a 

particular value or the calculations will become unstable. The first 

restriction Is that the fluid cannot flow through more than one cell in one 

single time step. This is because the finite difference expressions assume 
mass or momentum fluxes only between adjacent cells. Hence the time step At 
should be less than the minimum time taken for the fluid to pass through one 
cell, taken over all the cells In the grid. 


Therefore, 


At<min 


r 4r i Lz j l 

iK'./TVjlJ 


The At used for computation is usually taken to be 0.25-0.33 times that 
obtained from above [161. 

The second restriction arises from the fact that the fluid should not 
diffuse through more than one cell in a single time step. The expression 
obtained after performing a linear stability analysis [17] is 


fAt 



i>At 



< 


1 

2 
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7,7 List of Variables 

Symbol 

CONDHEIGHT 

thickr# 

zcond# 

CONVCT 

uLT#, uRT#, 
vFT#, vAT# 

tLwig#, tRwigtt, 
tFwig#, tAwig# 

CONVCU 

uLun#, uRun#, 
vFun#, vAun# 

uLwig#, uRwig# f 
uFwig#, uAwIg# 

CONVCV 

uLvn#, uRvn#, 
vFvn#, vAvn# 


or 


At < 


Zv 1 


Ar, 


Az 


Used In The Code 


Description 


Subroutine used to calculate condensate heights 
Condensate height at radial position radpos# 

Condensate height at location of cell center re# 

Subroutine which evaluates the convection term of 

the temperature equation; ^- ru ^ + 

r or oz 

Velocities uf 1 , u£ v" and v” _ at the left, 

L|T K* I i i 1 A f T 

right, front, and aft faces of the u-cell 

respectively 

Upstream values of temperature T^, ? R , Tp and T A at 

the left, right, front, and aft faces of the T-cell 

respectively 

Subroutine which evaluates the convection term of 

the r-momentum equation: — ®L ruu ^ + ££ vu ^ 

r dr dz 

Convecting velocities u” , u” , v™ and v? at 

L,u R,u F,u A,u 

the left, right, front, and aft faces of the u-cell 

respectively 

Convected velocities u^, u R , Up, and u^ at the left, 
right, front, and aft faces of the u-cell 

respectively 

Subroutine which evaluates the convection term of 

. .. 1 diruv ) 3(vv) 

the z-momentum equation: — g— + g- 

Convectlng velocities u£ y , u R y , Vp and at 

the left, right, front, and aft faces of the v-cell 

respectively 
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vLwigtt, vRwig#, 
vFwig#, vAwig# 

delt# 

DIFFNT 

alphaL#, alphaR#, 
alphaF#, alphaA# 

dtrL#, dtrR#, 
dtzF#, dtzA# 


DIFFNU 

nuLun#, nuRuntt, 
nuFun#, nuAun# 

durL#, durR#, 
duzF#, duzA# 


DIFFNV 

nuLvn#, nuRvn#, 
nuFvn#, nuAvn# 

dvrL#, dvrR#, 
dvzF#, d.vzA# 


Convected velocities v^, Vp, and v at the left, 
right, front, and aft faces of the v-cell 
respectively 


Time step, At 


Subroutine which computes the diffusion term of the 


a(«|I) 


. 1 d[ar^~) 

temperature equation: — gp dr + ^ dz 

Thermal diffusivities a?-, a" _, a” _, and a* at 

L| i f i • i i A y 1 

the left, right, front, and aft faces of the T-cell 
respectively 

_ . dT , dT. dT. . dT. ^ 

Derivatives W \ L , W \ R , ^| F and ^| A at the left, 


right, front 
respectively 


and aft faces of the T-cell 


Subroutine which computes the diffusion term of 
the r-momentum equation: p 

Kinematic viscosities i/ 1 , , vl and v? at 

L,u’ R,u* F,u A,u 

the left, right, front and aft faces of the u-cell 
respectively 

Derivatives g£| L , §£| R . |^| F and §£| A at the left, 
right, front and aft faces of the u-cell 
respectively 


Subroutine which computes the diffusion term of the 

, ,, i a(vrfi) alnfc 

z-momentum equation: — gp dr + dz 


Kinematic viscosities , iv and 

L,v R,v F,v 

v at the left, right, front and aft 

A, V 

faces of the v-cell respectively 

, .. dv . dv. 8v. . 8v, 

Derivatives jp| R . jj| f and jj| a at 

the left, right, front and aft faces of 

the v-cell respectively 
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epst#,epsv# 

errt# 

HEATFLUX 

hfluxr# 

qout# 

iflag7. 


INITIAL 

INPUTS 

cph2o#, 

rhoh2o# 

nuairtt, 

INTER PI 

tempr# 

tw# 

INTERP2 


OUTPUTS 


Convergence criteria for the pressure iteration 
procedure and temperature field respectively 

Error, (temp-tempn), after each time step 

Subroutine calculating heat fluxes 
Heat flux at radial position radpos# 

Heat flux at cell center location rc# 

Parameter used to determine whether outflow 
boundary conditions for velocity are to be 
applied or not. iflag7 ■ 0 for pressure 

iteration when outflow conditions are not 
applied, iflag7 = 1 otherwise 

Subroutine used to set the initial conditions 
for velocities and temperature 

Subroutine which reads model parameters from 
files "thermal. par" and "conden.par” 
kh2o#, nuh2o#, Specific heat at constant pressure, C , thermal 

, cpair#, kair#, conductivity, k f kinematic viscosity, v t and 

rhoair# density, p of water and moist air respectively 

Subroutine used to interpolate temperatures 

from radial positions radpos# to locations at 
cell centers re# 

Temperature at position radpos# 

Temperature at position rc# 

Subroutine used to interpolate condensate 
heights and heat fluxes from cell center 
locations rc# to specified radial positions 
radpos# 

Subroutine which writes out results to file 

"conden.par" 


52 



PROPERTY Subroutine defining physical properties at grid 

points 

cpll, k#, rho#, vise# Specific heat at constant pressure, thermal 

conductivity, density and kinematic viscosity 


p#, pn# 


Pressure p at present and previous time levels 
respectively 


PRSITR 

beta# 

delp# 

dij# 

omega# 


Subroutine which adjusts cell pressures and 
velocities until continuity is satisfied 
Geometric factor 

pressure adjustment A p 

Left hand side of the mass conservation 
equation 

relaxation factor w 


REGION 

delr#, delrb2#, delu#, 
re#, ru#, delz# f 
delzb2#, delv#, zc#, zv# 


Subroutine which calculates geometric 
parameters associated with the mesh 
The distances Ar, Ar/2, Au, rc, ru, Az, Az/2, 
Av, zc and zv respectively 


slope# 


T 


ij 


Term used to add the effects of convection, 
diffusion and source terms e.g. 


r u * “ 


slope where slope. 


,J_ 8j_ruT) ££v7).| 
- Pi-7 M 


r dr 


i a i 37. a# 37 . n 
1 8 ( ar 3—) a(a-a-), 

- — dr * dz + 


r dr 


pc. 


Dp 

Dt 


n 


SOURCT 


dpt#, udpr#, vdpz 
pL#, pR#, pF#, pA# 


uT#, vT# 


Subroutine which evaluates the source term of 
the temperature equation: 

P P D 

§f. u §7 and respectively 

Pressure values at the left, right, front and 
aft faces of the 7-cell respectively 
Velocities at the center of the 7-cell 


SOURCU 


Subroutine which evaluates the source term of 



dpr#, dvrF#, dvrA#, 


rhou# 


srcul#, srcu2#, srcu3# 


the r-momentum equation: 


I 3(vr^) 
r dr 8 


d[v~) 

Si 3r 



i ££ 

p dr 


- £ §7 and derivatives §pl F and If I a at the 

front and aft faces of the u-cell respectively 
Density at the center of the u-cell 

I and - 2i>-^ respectively 

r dr az ^2 r 


SOURCV 


dpz#, duzL#, duzR#, 


rhov# 

srcvl#, srcv2# 


Subroutine which evaluates the source term of 

„ I 8 ( vrjr~) . a(v|^) 
the z-momentum equation: — dz + ^ az 

1 dp 
p dz 

■ ^ If and derivatives Sii L and lfi R at the 

left and right faces of the v-cell respectively 
Density at the center of the v-cell 

■ «/ flu. a , 3v. 

I -^- vr dz and -a- V dz respectively 

r dr dz 


tamb#, tdewpt# 


Ambient and dew point temperatures of the moist 
air respectively 


temp#, tempn# 


Temperature T at the present and previous time 
levels respectively 


TBOUNDS 


Subroutine used to set the boundary conditions 
for temperature 


u#, un#, v#, vn# 


Velocities u and v at the present and previous 
time levels respectively 


vin# 

VBOUNDS 


Inlet velocity of moist air 

Subroutine which sets the velocity boundary 
conditions 
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7.8 Test Cases 

To check the condensate model three test cases were run, each with a 
different pressure wall temperature profile. For all the test cases, the 
following values were used for the model parameters: 

Ambient temperature = 294 K ■ 21 °C 

Dew point temperature « 288 K ■ 15 °C 

Coefficient of heat transfer for moist air » 5 W/(m Z K) 

Radius of pressure wall surface = 1 m 
inlet velocity = 0.005 m/s a 0.015 ft/s 

Physical properties of water were prescribed at the dew point temperature 
while those of moist air were prescribed at the ambient temperature. The 
temperature and condensate height distribution along the radial position for 
test cases I, II and III are shown in Figures 7.4, 7.5 and 7.6, and listed in 
Tables 7.2, 7.3, and 7.4, respectively. The results obtained appear to be 

physically reasonable. The authors could locate no appropriate experimental 
data In the literature for comparison with the calculated results. 
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I BUMPER 



100 mm 


MULTILAYER INSULATION (MLI) 

30 LAYERS DOUBLE ALUMINIZED 
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Fig. 1.1 Schematic drawing of impact specimen. 


ENTER DATA 


Test ID: 



Data Source: 


Test Date: 



Bunper Mat r l; 
Prs Wall Ilat'l 
Proj DlaatUr: 


6061-T6 


Bunper Thk: 


2219-T8 


Prs Wall Thk: 
lepact Angle: 


0.125 


Bunper Std-Off: 
Proj Hat' 1 : 
Proj Vel : 


1100 


Bunper Hole Kajor Axis: 
MLI Hole Diameter: 


Bunper Hole Minor Axis: 
Pressure Mall Hole Dianeter: 


< Add to Database > < Cancel this Data Entry > < Exit Progran > 


Fig. 2.3.1 Data entry window for adding records to the impact data file. 
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Number of Data Record To Remove: | 



< OK to Removs > < Quit > 


Fig. 2.3.2 Data entry window for deleting records from the Impact data 

file. 


VIEW DATA 


Test ID: 


1012 

Data Source: 

HSFC 

Test Date: 

05/08/9 


Bumper Hat'l 
Prs Wall Hat'l: 
Proj Diameter: 


6061-T6 


Bumper Thk: 


.08 


2219-T8 


Pra Wall Thk: 


.125 


Bumper Std-Off: 
Proj Hat'l: 


1100 


.313 

Impact Angle: 

0 

Proj Vel: 



Bumper Hole Major Axis 
KLI Hole Diameter: 


.729 


Bumper Hole Minor Axis: 


.729 


2.2 


Pressure Wall Hole Diameter: 


.375 


< Heart Data Record > 


< Exit Program > 


Fig. 2.3.3 Data entry window for viewing records in the impact data file. 


EDIT IMPACT PAR 


Impact Data File: 


mli.dat 


Material Data File: 


material.dat 


Impact Results File: 


impact. out 


Bumper Hat'l 
Prs Wall Hat'l: 
Proj Diameter 


6061-T6 


Bumper Thk: 


0.08 


2219-T8 


Prs Wall Thk: 


0.125 


Bumper Std-Off: 
Proj Hat'l: 


1100 


0.25 

Impact Angle: 

45 

Proj Vel: 

7.1 


< Save Changes and Exit > 


< Exit Program > 


Fig. 2.3.4 Data entry window for editing the impact parameters file. 
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IMPACT DAMAGE 1 

Bumper Hole Major Diameter: 

Bumper Hole Minor Diameter: 

MLI Hole Diameter: 

Pressure Wall Hole Diameter: 

< OK > 

Fig. 2.3.6 Data window for displaying impact results. 



1 test 

CENTERLINE 

j BUMPER SIDE 

1 

1 

1 

1 

1 

BUMPER DATA 

Centerline Temp: 
115.9 

Border Tenp : 
124.1 

Hole Dianeter: 
0.01643 

1 

* 

MLI Hole Dian: 
0.06977 

1 

1 

PRESS . WALL DATA 

Center 1 ine Tenp : 
66.0 

Border Tenp : 

66.3 

1 

1 

1 

1 

1 

1 

1 PRESSURE HALL SIDE 

I 

Press <SPACE BAR> to Continue 

MAX CONDEN THICK 
0.00 


Fig. 2.3.7 Graphics screen showing thermal and condensation calculation 
results. 
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LEGEND 

POINTS WHERE THE 
■ DEPENDENT VARIABLE 
HAS BEEN MEASURED 

POINT WHERE A 
^PREDICTION IS REQUIRED 
(TARGET POINT) 

• INTERPOLATION POINTS 
AMIN/MAX POINTS 


INDEPENDENT VARIABLE 1 


Fig. 3.1 Technique for selecting interpolation point locations for the 
case of two independent variables. 


INDEPENDENT 
VARIABLE 2 



LEGEND 

„ MEASURED VALUES OF 
DEPENDENT VARIABLE 


POINT WHERE A 
• PREDICTION IS REQ'D 
(TARGET POINT) 


INDEPENDENT VARIABLE I* 


Fig. 3.2 Interpolation scheme for equally spaced data points. 


INDEPENDENT 
VARIABLE 2 



LEGEND 

MEASURED VALUES OF 
DEPENDENT VARIABLE 

POINT WHERE A 
| PREDICTION IS REQ'D 
(TARGET POINT) 


INDEPENDENT VARIABLE I 


Fig. 3.3 Interpolation scheme for unequally spaced data points. 
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BUMPER SIDE 




N+1 th LAYER 



OUT OF PLANE 
RADIATION AND 
CONDUCTION 


N th LAYER 


IN-PLANE 
CONDUCTION 


N-1 th LAYER 


Fig. 6.3 Schematic drawing of heat flow into and out of a typical node. 
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Fig. 6.5 Ring to ring view factor geometry. 



MLI HOLE 



Fig. 6.6 Heat flux from rings of the first MLI layer to a node in a pressure 
wall ring. 


MLI HOLE 


HEAT FLUX FROM SPACE 



Fig. 6.7 Schematic drawing Illustrating the method of calculating heat flux 
to the pressure wall from the bumper and space environment through the MLI 
hole. 
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Jmix 


jmavl 



l-l 2 


FRONT 


imix-1 Imtx 


RIGHT 


Fig* 7.2 Staggered mesh system showing locations of velocities and 
temperatures. 
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Fig. 7.3 Dimensions of different cells used in the computation of u,v and t. 










Test Run 


Calculation Results Data Pile Name: thermal. out 

Initial Values File Name: thermal.ini 


MLI hole diameter: 

MLI stand off: 

Estimated pressure wall temperature: 
Estimated bumper temperature: 

Temperature conversion factor one: 
Temperature conversion factor two: 

Number of MLI layers: 

Radius of area modeled: 

Pressure wall thickness: 

MLI layer thickness: 

Beta cloth thickness: 

Bumper thickness: 

Bumper stand off: 

Space Thermal Radiation Plux: 

Thermal conductivity of pressure wall: 
Thermal conductivity of MLI: 

Heat transfer coefficient of Dacron Netting; 
Thermal conductivity of beta cloth: 

Thermal conductivity of the bumper: 
Emissivity of pressure wall: 

Emissivity of MLI: 

Emissivity of beta cloth: 

Emissivity of outer surface of bumper: 
Emissivity of inner surface of bumper: 
Stefan-Boltzmann constant: 

Maximum number of iterations for eech mesh: 
Convergence Factor: 

Initial Number of Nodes: 

Maximum Number of Nodes: 

Bumper Hole Diameter 
Module Air Temperature: 

Module air dew point temperature: 

Convective heat transfer coefficient: 
Condensate density: 

Condensate kinematic viscosity: 

Condensate thermal conductivity: 

Condensate constant prsssure specific heat: 
Module air density: 

Module air kinematic viscosity: 

Module air thermal conductivity: 

Module sir constant pressure specific heat: 


.069767 

.0508 

295 

100 

-459.67 

1.8 

20 

.5 

.003175 

.00000635 

.0000508 

.001524 

.1016 

431 

130 

50 

1.0687 

5 

115 

.06 

.06 

.94 

.94 

.14 

.000000056697 

10000 

.0001 

10 

10 

.016427 

295 

290 

5 

1000.52 

.000001006 

.597 

4181.8 

1.1774 

.00001568 

.02624 

1005.7 


Calculations Stopped by User Before Convergence! 


Pinal Nodal 

Temperatures : 


Node No. 

Pressure Mall 

Bumper 

1 

6 . 603D+01 

1 . 159D+02 

2 

6.603D+01 

1 , 159D+02 

3 

6 • 603D+01 

1 • 159D+02 

4 

6 . 609D+01 

1 . 168D+02 

5 

6 ♦ 615D+01 

1.182D+02 

6 

6 . 621D+01 

1 . 198D+02 

7 

6 . 6260+01 

1.216D+02 

8 

6 . 629D+01 

1.231D+02 

9 

6 .631D+01 

1 . 241D+02 

10 

6 . 63 1D+01 

1 . 241D+02 

Node No. 

Condensate Thickness: 

1 

0.000D+00 


2 

0 . 00OD+00 


3 

0 . OOOD+OO 


4 

0. 000D+00 


5 

0. 000D+00 


6 

0. 000D+00 


7 

0. 000D+00 


8 

0. 000D+00 


9 

0. 000D+00 


10 

0. 000D+00 



Table 2.3.1 Typical thermal and condensation calculations results file. 
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Order In Which 
Daia Seu 
Axe Tested For 
Linear Independence 

Data Points Selected To Form Set Of Four 
(Ordered From Closest To Prediction Point To Farthest) 

1 2 

3 

4 

5 

6 

7 

I 

1 

2 

3 

4 




2 

1 

2 

3 


4 



3 

1 

2 

3 



4 


4 

1 

2 

3 




4 | 

5 

1 

2 


3 

4 



6 

1 

2 


3 


4 


7 

I 

2 


3 



4 | 

8 

1 

2 



3 

4 


9 

1 

2 



3 


4 

10 

1 

2 




3 

4 

11 

1 


2 

3 

4 



12 

I 


2 

3 


4 


13 

1 


'S 

3 



4 | 

14 

1 


2 


3 

4 


13 

I 


2 


3 


4 

16 

1 


2 



3 

4 

17 

I 



2 

3 

4 


18 

l 



2 

3 


4 

19 

l 



2 


3 

4 

20 

1 




2 

3 

4 


Table 4.1 Scheme for Selecting Four Data Point Sets from the Closest Seven 
Nodes for Damage Function Coefficient Determination. 


Equation 

* 

r * 

% 

Continuity 

1 

0 

0 

r-Momentum 

u 


TF + rVr 8r * Si ar 

2-Momentum 

V 

I* 

-5i*? iH*’ • S* 

Energy 

h 

k/C 

E- 

§f 


Table 7.1 Summary of Equations. 



Radial Position 
(mi 

Temperature 

(K) 

Heat Flux 
(W/m*J 

Condensate 
Height (ml 

o.ooooo 

283.00000 

167.22690 

0.01851 

0.05000 

283.30000 

160.47165 

0,01749 

0.10000 

283.70000 

153.71641 

0.01646 

0.15000 

284.20000 

146.96116 

0.01544 


0.20000 

284.80000 

137.19782 

0.01357 



0.25000 

285.50000 

127.43448 

0.01171 



0.30000 

286.50000 

111.61815 

0.00710 

Table 7.2 

Results of condensate test case Li 

0.35000 

287.60000 

95.80183 

0.00249 


0.40000 

288.90000 

57.40092 

0.00000 



0.45000 

290.20000 

19.00000 

0.00000 


i 

0.50000 

291.40000 

14.50000 

0.00000 


0.55000 

292.00000 

10.00000 

0.00000 



0.60000 

292.20000 

9.00000 

0.00000 



0.65000 

292.40000 

8.00000 

0.00000 



0.70000 

292.70000 

6.25000 

0.00000 



0.75000 

293.10000 

4.50000 

0.00000 



0.80000 

293.50000 

3.00000 

0.00000 



0.85000 

293.70000 

1.50000 

0.00000 


- 

0.90000 

293.80000 

1.00000 

0.00000 


- 

0.95000 

293.90000 

0.50000 

0.00000 



1.00000 

294.00000 

0.25000 

0.00000 




Radial Position 
(m) 

Tempearture 

<K1 

Heat Flux 
(W/m 1 ) 

Condensate 
Height (ml 

0.00000 

286.00000 

122.26894 

0.01084 

0.05000 

286.30000 

115.49804 

0.00879 

0.10000 

286.70000 

108.72715 

0.00674 

0.15000 

287.20000 

101.95625 

0.00468 

0.20000 

287.80000 

64.72812 

0.00234 

0.25000 

288.50000 

27.50000 

0.00000 — . . _ _ 

0.30000 

289.50000 

22.25000 

0.00000 laOle /.o 

0.35000 

290.60000 

17.00000 

0.00000 

0.40000 

291.90000 

13 00000 

0.00000 

0.45000 

292.20000 

9.00000 

0.00000 

0.50000 

292.40000 

7.00000 

0.00000 

0.55000 

293.00000 

5.00000 

0.00000 

0.60000 

293.10000 

4.50000 

0.00000 

0.65000 

293.20000 

4.00000 

0.00000 

0.70000 

293.30000 

3.37500 

0.00000 

0.75000 

293.45000 

2.75000 

0.00000 

0.80000 

293.55000 

2.25000 

0.00000 

0.85000 

293.65000 

1.75000 

0.00000 

0 90000 

293.80000 

1.12500 

0.00000 

0.95000 

293.90000 

0.50000 

0.00000 

1.00000 

294.00000 

0.25000 

0.00000 


Results of condensate test case II: 


Radial Position 
(mi 

Temperature 

(K) 

Heat Flux 
(W/m*) 

Condensate 
Height (m) 

0.00000 

288.10000 

29.50000 

0.00000 

0.05000 

288.30000 

28.50000 

0.00000 

0.10000 

288.70000 

26.50000 

0.00000 

0.15000 

289.20000 

24.00000 

0.00000 

0 20000 

289.80000 

21.00000 

0.00000 

0.25000 

290.50000 

17.50000 

0.00000 Table 7.4 

0.30000 

291.50000 

12.50000 

0.00000 

0.35000 

292.60000 

7.00000 

0.00000 

0.40000 

292.90000 

5.50000 

0.00000 

0.45000 

293.20000 

4.00000 

0.00000 

0.50000 

293.40000 

3.00000 

0.00000 

0.55000 

293.55000 

2.25000 

0.00000 

0.60000 

293.65500 

1.72500 

0.00000 

0.65000 

293.75900 

1.20500 

0.00000 

0.70000 

293.88000 

0.60000 

0.00000 

0.75000 

293.91000 

0.45000 

0.00000 

0.80000 

293.93000 

0.35000 

0.00000 

0.85000 

293.95500 

0.22500 

0.00000 

0.90000 

293.97500 

0.12500 

0.00000 

0.95000 

293.99000 

0.05000 

0.00000 

1.00000 

294.00000 

0.00000 

0.00000 


Results of condensate test case III 
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